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Kurzfassung 


Die Funktionserfüllung vieler technischer Systeme lässt sich auf den gleitenden Kontakt 
zweier oder mehrerer Körper zurückführen. Hier können zum Beispiel Bremsen, 
Kupplungen und Gleitlager genannt werden, die praktisch in jedem mobilen oder 
stationären Antriebssystem zu finden sind. Für die Weiterentwicklung dieser reibungs- 
behafteten Systeme - sei es aus funktionaler oder ökonomischer Sicht - ist entsprechend 
ein möglichst umfassendes Verständnis der zugrundeliegenden Zusammenhänge im 
Kontakt notwendig. 

Die vorliegende Arbeit beschäftigt sich daher mit der Fragestellung, inwieweit der 
Reibwert für den Kontakt zweier gleitender metallischer Körper in Abhängigkeit 
verschiedener physikalischer Parameter bei trockener Reibung berechnet werden kann. 
Darüber hinausgehende Untersuchungen zu den Auswirkungen des berechneten 
Reibwerts im Kontext von reibungserregten Schwingungen schließen sich an. 

Dafür wird ein thermomechanisches Modell für den Kontakt zweier gleitender Kör- 
per mit rauen Oberflächen entwickelt und für verschiedene Kontaktkonfigurationen 
ausgewertet. Die Modellbildung basiert dabei auf der Halbraumannahme und den 
quasi-statischen Bilanzgleichungen für einen thermoelastischen Festkörper. Das Simu- 
lationsmodell ermöglicht die Berechnung verschiedener kontaktmechanischer Größen 
wie beispielsweise den Verschiebungen, den Spannungen, dem Temperaturfeld, der 
tatsächlichen Kontaktfläche und insbesondere dem stochastischen Reibwert. Es ist 
außerdem in der Lage, die Einflüsse von temperaturabhängigen Materialparametern 
abzubilden. Die vorliegende Arbeit schließt damit eine in der Literatur noch vorhandene 
Modellierungslücke und kann als Erweiterung der bisher existierenden Kontaktmodelle 
angesehen werden. Weiterführende Untersuchungen hinsichtlich den zugrundeliegen- 
den Instabilitätsmechanismen bei reibungserregten Schwingungen unter Einbezug der 
Erkenntnisse aus den thermomechanischen Simulationen schließen sich an. Hierfür 
wird der berechnete Reibwert zunächst geeignet parametrisiert und anschließend in 
ein Minimalmodell eingearbeitet. Die Stabilitätsanalysen zeigen bisher noch nicht 
in der Literatur bekannte stabilisierende und destabilisierende Auswirkungen des 
stochastischen Reibwerts auf Systeme mit tribologischen Kontakten. 
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1 Einleitung 


Selten lassen sich komplexe physikalische Zusammenhänge durch triviale mathemati- 
sche Ausdrücke beschreiben. Dass dies bei der Reibkraft 


F, DÉI (1.1) 


so gut gelingt, ist deshalb umso erstaunlicher. Die Reibkraft F, beschreibt den Wider- 
stand, der einer Relativbewegung von zwei sich berührenden Körpern entgegenwirkt 
und als Produkt des Reibwerts u mit der Normalkraft F„ berechnet werden kann. Die 
Modellierung der Reibkraft durch Gleichung (1.1) ist heutzutage in jedem kommerzi- 
ellen Simulationsprogramm hinterlegt und gemeinhin als Coulomb’sches Reibgesetz 
bekannt. Allerdings gelingt es bisher trotz intensiver Bemühungen nicht, den Reibwert u 
für beliebige Kontaktkonfigurationen zuverlässig vorherzusagen. Aus diesem Grund 
beschäftigt sich auch die vorliegende Arbeit mit dieser Thematik. Sie reiht sich damit 
in die Chronologie tribologischer Forschungsarbeiten ein, deren Ursprünge bis in 
die Renaissance zurückreichen. Der Fokus der Arbeit liegt dabei auf dem gleiten- 
den Kontakt metallischer Körper mit trockener Reibung, wie sie im Maschinenbau 
häufig Anwendung finden sowie den damit verbundenen Auswirkungen hinsichtlich 
reibungserregter Schwingungen. Warum die Berechnung des Reibwerts u nicht nur aus 
akademischer Sicht von Interesse ist, wird in Unterkapitel 1.1 motiviert. Anschließend 
folgt ein Literaturüberblick in Unterkapitel 1.2, der sich größtenteils auf den gleitenden 
Kontakt zweier Körper mit trockener Reibung konzentriert. Kapitel 1 wird nach einer 
Synthese des Forschungsstandes durch eine Beschreibung von der Struktur und der 
Zielsetzung der Arbeit in Unterkapitel 1.3 abgeschlossen. 


1.1 Motivation 


Die Funktionserfüllung vieler technischer Systeme lässt sich auf den gleitenden Kontakt 
zweier oder mehrerer Körper zurückführen. Hier können zum Beispiel Bremsen, 
Kupplungen und Gleitlager genannt werden, die sich im Allgemeinen in jedem mobilen 
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oder stationären Antriebssystem finden lassen. Eine Modellierung der Reibkräfte in 
den jeweiligen tribologischen Kontakten durch Gleichung (1.1) ist in Anbetracht der 
dort wirkenden komplexen physikalischen Prozesse vorteilhaft. Entsprechend groß ist 
jedoch auch die Anzahl der Parameter, die einen Einfluss auf den Reibwert oder den 
zugrundeliegenden tribologischen Kontakt haben [241]. Experimente und Simulationen 
zeigen unter anderem eine Abhängigkeit des Reibwerts von den beteiligten Werkstoffen, 
der Gleitgeschwindigkeit, der Temperatur und der Oberflächenstruktur. 

Aus der Perspektive des Entwicklungsingenieurs lässt sich davon neben der Werkstoff- 
auswahl und der Betriebstemperatur meist nur die Struktur der Kontaktflächen gezielt 
beeinflussen. Durch moderne Fertigungsverfahren sind mittlerweile Bearbeitungen der 
Oberflächen bis in den nm-Bereich möglich [12, 90], eine Gestaltabweichung in Form 
von Oberflächenrauigkeiten ist dabei dennoch immer vorhanden. Veränderungen der 
Oberflächenstruktur können zu einer deutlichen Reduktion des Reibwerts führen [90, 
221], die Tragfähigkeit erhöhen [37] oder die Dynamik des übergeordneten technischen 
Systems verbessern [45, 100, 109, 126, 249]. Es lässt sich aber ebenso das Verschleiß- 
verhalten beeinflussen [61, 90], der Transport von Schmiermittel steuern [127] oder 
der Wärmeübertrag und die elektrische Leitfähigkeit zweier sich berührender Körper 
maßgeblich steigern [61]. In allen Fällen spielt die tatsächliche Kontaktfläche zwischen 
den beiden Kontaktkörpern eine wesentliche Rolle. 

Neben den beschriebenen Einflüssen auf die Funktionalität von technischen Systemen 
lassen sich auch ökonomische Aspekte der Reibkraft feststellen. Für die Energieverluste 
aufgrund von unerwünschter Reibung müssen jährlich hunderte Millionen Euro 
aufgewendet werden [240]. Studien zu den Reibungsverlusten in Personenkraftfahr- 
zeugen [76, 103], in Lastkraftwagen und Bussen [104] sowie in der Bergbauindus- 
trie [105] zeigen, dass Reibungsverluste signifikant zum Energieverbrauch beitragen. 
Schätzungsweise 20 % des weltweiten Energiebedarfs werden aufgrund von parasitären 
Reibungsverlusten verursacht [106]. Diese tragen somit erheblich zu den globalen 
Treibhausgasemissionen bei. 

Für die Weiterentwicklung von technischen Systemen mit Reibung - sei es aus funktiona- 
ler oder ökonomischer Sicht - ist in allen Fällen ein möglichst umfassendes Verständnis 
der zugrundeliegenden Zusammenhänge nötig. Dass dadurch grundsätzlich ein großes 
Verbesserungspotential vorhanden ist, zeigen die aufgeführten Veröffentlichungen. 


1.2 Stand der Forschung - Von der Renaissance bis in 
das 21. Jahrhundert 


Die folgenden Ausführungen konzentrieren sich auf den trockenen Kontakt rauer Kör- 
per und orientieren sich für die historischen Entwicklungen zwischen der Renaissance 
und der Mitte des 20. Jahrhunderts an den Zusammenfassungen [24, 68, 77, 162, 208, 
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250]. Für den Zeitraum ab dem 19. Jahrhundert bis heute fließen zunehmend auch 
weitere, als relevant erachtete Veröffentlichungen ein, sodass sich schlussendlich ein 
möglichst umfassendes Gesamtbild des aktuellen Forschungsstandes ergibt. 


Die ersten wissenschaftlichen Untersuchungen auf dem Gebiet der Tribologie werden 
in die Zeit der Renaissance datiert und sind Leonardo da Vinci (*1452 — +1519) zuge- 
schrieben. In seinen Experimenten stellt da Vinci eine Proportionalität des Widerstands 
gegen den Beginn des Gleitens zum Gewicht des Körpers fest, findet jedoch keinen 
Einfluss der scheinbaren Kontaktfläche. Er versteht dennoch die Oberflächenstruktur 
als Ursache des Widerstands und schlussfolgert eine kleiner werdende Reibkraft mit 
glatter werdender Oberfläche. 

Erst mit dem Beginn des 18. Jahrhunderts werden weitere Untersuchungen von 
Guillaume Amontons (*1663 — +1705) durchgeführt. Er formuliert unabhängig von 
da Vinci Gesetzmäßigkeiten der Reibung, die aber im Einklang mit dessen Beobach- 
tungen stehen. Zudem sieht Amontons ineinandergreifende Asperiten und eine damit 
notwendige Hebungsarbeit als Ursache für den Widerstand gegen den Beginn des 
Gleitens an und rückt damit die Interaktion der Oberflächenstrukturen als Ursprung 
für Reibung in den Vordergrund. Eine Vorstellung, die ebenfalls von Philippe de la Hire 
(*1640 - +1718) geteilt wird, der mit seiner Arbeit die Resultate Amontons bestätigt. 
Insbesondere vermutet de la Hire, dass für den Beginn des Gleitens unter anderem 
auch die Zerstörung einzelner Asperiten notwendig sein könnte. Aus heutiger Sicht 
kann dies als plastische Verformung und damit permanente Oberflächenveränderung 
interpretiert werden. Mit John Theophilus Desaguliers (*1683 - +1744) entsteht die neben 
den ineinandergreifenden Asperiten konzeptionell neue Vorstellung der Adhäsion als 
Ursprung für die Reibkraft. Durch ihn entwickelt sich darüber hinaus erstmals der 
Gedanke, dass die Reibkraft bei immer glatter werdenden Oberflächen wieder steigen 
statt fallen könnte. Frühe Modelle von rauen Oberflächen entstehen. Bernard Forrest de 
Belidor (*1697 - +1761) repräsentiert diese als Ansammlung kugelförmiger Asperiten, 
während Leonhard Euler (*1707 - +1783) eine Modellierung der Asperiten in der Form 
von Dreiecken wählt. Euler prägt zudem die Nutzung des griechischen Buchstabens u 
für den Reibwert und unterscheidet zwischen Haft- und Gleitreibung. Mit den 
umfangreichen Arbeiten von Charles Augustin Coulomb (*1736 - +1806) wird schließlich 
ein neues Kapitel in der Tribologie aufgeschlagen und die Forschungsbemühungen 
erhalten ein wesentlich komplexeres und größeres Ausmaß. Coulomb untersucht unter 
anderem den Einfluss von der zeitlichen Dauer der Haftphasen sowie von verschiedenen 
Umgebungseinflüssen wie beispielsweise Temperatur und Luftfeuchtigkeit. Außerdem 
interpretiert Coulomb die Asperiten als verformbare Borsten. Erste Versuche, die 
Abhängigkeiten des Reibwerts in analytischen Zusammenhängen auszudrücken, gehen 
ebenfalls auf ihn zurück. 
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Damit besteht bereits zu Anfang der industriellen Revolution Ende des 18. Jahrhunderts 
eine recht umfangreiche Vorstellung für die Ursachen und Zusammenhänge der 
Reibkraft, von denen viele Einzug in die Forschungen jüngerer Zeit gefunden haben. 
Im 19. Jahrhundert sind die Entwicklungen in der Tribologie größtenteils durch die 
Erforschung der Fluidmechanik geprägt. Mit der nach Osborne Reynolds (*1842 - +1912) 
benannten Reynoldsgleichung und den Experimenten von Richard Stribeck (*1861 - 
+1950) werden im Bereich der geschmierten Kontakte große Fortschritte erzielt. Der als 
Stribeck-Kurve bekannte Zusammenhang zwischen Reibkraft und Gleitgeschwindigkeit 
bei Gleitlagern wird mittlerweile auch oft im Kontext von geschwindigkeitsabhängigen 
Reibkräften in trockenen Kontakten verwendet. Untersuchungen zu trockenen Kon- 
takten fokussieren sich in dieser Zeit hauptsächlich auf den Übergang der Zustände 
Haften und Gleiten. Aus zahlreichen Experimenten schließen Henry Charles Fleeming 
Jenkin (*1833 — +1885) und James Alfred Ewing (*1855 - +1935) schließlich, dass der 
Ubergang des einen Zustandes in den anderen kontinuierlich ist. Als Ursache fiir die 
Reibkraft wird allerdings immer noch die fiir das Ubereinandergleiten der Asperiten 
notwendige Hebungsarbeit angesehen. 

Das 19. Jahrhundert beinhaltet ebenso einige Veröffentlichungen, denen aus heutiger 
Sicht eine große Bedeutung zukommt. Hier sind beispielsweise Joseph Valentin Boussi- 
nesq (*1842 — +1929) und Valentino Cerutti (*1850 — +1909) zu nennen [119], denen die 
Herleitung der Fundamentallösungen für normale und tangentiale Einzelkräfte auf der 
Oberfläche eines Halbraums gelingt. Ihre Ergebnisse finden in einer Vielzahl von aktuel- 
len Publikationen Anwendung und sind ebenfalls Bestandteil der vorliegenden Arbeit. 
Aber auch die Ergebnisse von Heinrich Rudolf Hertz (*1857 - +1894) [92] zum elastischen 
Kontakt glatter, nicht-konformer Körper sind aufzuführen, die bis heute viele Forscher 
zu neuen Erkenntnissen inspirieren. Zum einen können Johnson-Kendall-Roberts [118] 
und Derjaguin-Muller-Ioporov [64, 173] durch die Ergänzung des Hertz-Kontaktes 
um Adhäsionskräfte deren Einfluss auf die Kontaktfläche zeigen. Zum anderen sei 
die, jeweils unabhängig von Cesare Cattaneo (*1912 — +1943) und Raymond David 
Mindlin (*1906 - +1987) [169, 170] durchgeführte, Erweiterung des Hertz-Kontaktes 
um tangentiale Belastungen erwähnt. Hierdurch wird die auch als partielles Gleiten 
bezeichnete Übergangsphase zwischen den Zuständen Haften und Gleiten besser 
verstanden und gleichzeitig ein Zusammenhang der Reibkraft zur Verschiebung und 
insbesondere zur Vorgeschichte des Kontaktes nachgewiesen. Als partielles Gleiten 
wird ein Zustand bezeichnet, bei dem Teilgebiete der Kontaktfläche haften, während 
andere Teilgebiete derselben Kontaktfläche noch oder schon gleiten. Die Ergebnisse 
von Cattaneo-Mindlin decken sich mit experimentellen Beobachtungen [117], stehen im 
Einklang mit den Schlussfolgerungen von Jenkin und Ewing und werden auch aktuell 
erneut diskutiert [75]. 

In der ersten Hälfte des 20. Jahrhunderts verändert sich durch die Arbeiten von 
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Frank Philip Bowden (*1903 — +1968) und David Tabor (*1913 - +2005) [36] schließlich 
dauerhaft die Sicht auf die Vorgänge trockener Reibung. Die tatsächliche, anstatt 
der scheinbaren Kontaktfläche und die bereits von Desaguliers angedachte Rolle der 
Adhäsion werden als zentrale Ursachen für den oft linearen Zusammenhang von 
Normal- und Reibkraft erkannt. Den Ausführungen von Bowden und Tabor folgend, 
plastifizieren die in Kontakt tretenden Asperiten aufgrund der kleinen tatsächlichen 
Kontaktfläche und den damit lokal sehr hohen Drücken sofort. Sie schließen daraus, 
dass die Belastung direkt proportional zur tatsächlichen Kontaktfläche ist. Bei ebenfalls 
angreifenden, tangentialen Belastungen scheren die Asperiten bei Überschreitung der 
Versagensscherspannung ab, womit die bereits von da Vinci und Amontons festgestellte, 
heutzutage als Coulomb’sches Reibgesetz bekannte Proportionalität zwischen Normal- 
und Reibkraft und ebenso die Unabhängigkeit von der scheinbaren Kontaktfläche für 
den rein plastischen Lastfall begründet werden kann. Bowden und Tabor schlussfolgern 
zudem, dass die Adhäsion den im Vergleich zu pflügenden Verlusten oder zur 
elastischen Hysterese dominierenden Anteil der Reibkraft bei Kontaktkörpern mit 
Härten in derselben Größenordnung und nominell glatten Oberflächen ausmacht. Zu 
diesem Ergebnis kommen auch neuere, numerische Untersuchungen [28]. 

Seit den Arbeiten von Bowden und Tabor lässt sich eine zunehmende Verzweigung der 
Forschungsgebiete in immer komplexer werdende Teilgebiete feststellen, für die eine 
chronologische Aufarbeitung nicht weiter zielführend scheint und die in ihrer Gesamt- 
heit auch nicht mehr zu erfassen ist. Für einen Großteil der Forschungsbemühungen im 
Kontext der Tribologie steht aber ab der zweiten Hälfte des 20. Jahrhunderts sowohl die 
tatsächliche Kontaktfläche als auch die Verteilung der einzelnen Kontaktgebiete im Vor- 
dergrund. Konnte für den rein plastischen Kontakt nun die Proportionalität zwischen 
Normal- und Reibkraft argumentiert werden, gelingt dies für den rein elastischen Kon- 
takt erst in der zweiten Hälfte des 20. Jahrhunderts. Durch fortschreitende Messtechnik 
sind nun Vermessungen rauer Oberflächen möglich. Ernest James Abbott und Floyd 
Firestone charakterisieren als erste die Höhenwerte durch die nach ihnen benannte 
Abbott-Firestone Kurve. Es findet zunehmend eine Interpretation der Oberflächen als 
zweidimensionales Zufallsfeld statt [152, 153] und für die Charakterisierung werden 
unter anderem Verteilungsdichte und Leistungsdichtespektrum herangezogen [178]. 
Durch diese Erkenntnisse modellieren Greenwood und Williamson [88] die Asperiten 
von rauen Oberflächen als Ansammlung sich nicht beeinflussender Halbkugeln, deren 
Höhenwerte im Gegensatz zu den Vorstellungen von de Belidor einer Normalverteilung 
folgen. Unter Verwendung der Ergebnisse von Hertz berechnen sie eine mit der 
Normalkraft nahezu proportional ansteigende, tatsächliche Kontaktfläche. Mit der 
Annahme einer zur Kontaktfläche proportionalen Versagensscherspannung kann damit 
eine Begründung für die als Coulomb’sches Reibgesetz bekannte Proportionalität 
zwischen Normal- und Reibkraft im rein elastischen Lastfall gefunden werden. Der 


1 Einleitung 


Modellierungsansatz von Greenwood und Williamson verursacht eine hohe Anzahl an 
Folgearbeiten und erfreut sich auch heute noch großer Beliebtheit. An dieser Stelle sind 
beispielsweise die Erweiterungen von Bush, Gibson und Thomas [41] sowie von Fuller 
und Tabor [81] zu nennen, die die ursprüngliche Modellierung von Greenwood und 
Williamson um repräsentative parabolische Asperiten und Adhäsionskräfte erweitern. 
Nach Aussagen des Autors Greenwood [89] ist der grundsätzliche Modellierungsge- 
danke jedoch fehlerhaft. Er verweist auf Arbeiten von John Frederick Archard [8], der 
bereits einige Jahre vorher raue Oberflächen als selbstähnliche Asperiten mehrerer 
Größenskalen modelliert und für den rein elastischen Fall eine näherungsweise 
Proportionalität zwischen tatsächlicher Kontaktfläche und Kontaktlast nachweist. In die 
Reihe der analytischen Arbeiten für den rein elastischen Kontakt rauer Körper gliedert 
sich aus jüngerer Zeit Bo Persson [63, 201] ein, dessen Modell als Eingangsgröße 
das Leistungsdichtespektrum der Oberflächenhöhen beinhaltet und ebenfalls auf 
eine Proportionalität zwischen tatsächlicher Kontaktfläche und Kontaktlast schließen 
lässt. Dabei liegen die berechneten Proportionalitätsfaktoren von Bush, Gibson und 
Thomas sowie Persson zwar in derselben Größenordnung, die genauen Abhängigkeiten 
werden aber immer noch diskutiert [42, 63, 107, 108, 188]. Zudem lässt sich bereits 
für nicht-normalverteilte Oberflächenhöhen keine Proportionalität mehr zwischen 
tatsächlicher Kontaktfläche und Kontaktlast nachweisen [52, 125, 131]. 

Sowohl die rein elastischen als auch die rein plastischen Modellierungsansätze sind 
im Rahmen ihrer Annahmen zu verstehen, reale Kontakte stellen im Allgemeinen eine 
Mischform der beiden Grenzfälle dar. Dennoch basieren die gewonnenen Erkenntnisse 
und Schlussfolgerungen allesamt auf Oberflächenrauigkeiten. Wohl nicht zuletzt aus 
diesem Grund gibt es eine Vielzahl an analytischen und durch die seit Ende des 20. 
Jahrhunderts breit verfügbare Rechenleistung auch numerischen Untersuchungen zu 
elastischen, plastischen und elastisch-plastischen Kontakten rauer Körper. Dabei hat 
sich in diesem Bereich mittlerweile neben einigen Simulationen mithilfe der Finite- 
Elemente-Methode [107, 108] die Rand-Element-Methode weitestgehend durchgesetzt. 
Es entstehen Arbeiten zu quasi-statischem Normalkontakt mit dem Einfluss von 
Adhäsion [174, 175], zur Modellierungen von Plastizität mit und ohne Verfestigung [53, 
91, 113, 156, 179, 180, 238], zu beschichteten Körpern [195, 196, 254, 255, 261], zu 
nicht-normalverteilten Oberflächenhöhen [52, 125, 131], zu partiellem Gleiten und 
tangentialer Belastung [38, 48, 82, 85, 98, 189, 190, 254, 255, 257], zu den Einflüssen 
der spektralen Komponenten der Oberflächenhöhen [256, 259, 260, 263] sowie zu 
Berechnungen der resultierenden Kontaktsteifigkeiten [43, 168] und Spannungen [145]. 
Auch das von Archard [7] eingeführte und nach ihm benannte Verschleiß-Gesetz 
sei erwähnt, dessen kontaktmechanischer Ursprung ebenfalls Gegenstand aktueller 
Forschungsbemühungen ist [80]. Die Mehrheit der aufgeführten Arbeiten kommt jedoch 
nicht mehr ohne numerische Hilfsmittel aus. Außerdem gehen alle Arbeiten entweder 
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von dem statischen Normalkontakt zweier rauer Körper aus, für den die zweiseitigen 
Rauigkeiten der Kontaktkörper zu einer resultierenden Rauigkeit zusammengefasst 
werden können, oder die Rauigkeit eines Kontaktkörpers wird vernachlässigt. Eine 
konsequente Modellierung im Sinne einer transienten äquivalenten rauen Oberfläche 
für gleitende Kontakte wurde nicht gefunden. Vielleicht gerade deswegen wird auf 
ein noch fehlendes Verständnis von gleitenden Kontakten hingewiesen [46]. Die erst 
kürzlich durchgeführte Contact Mechanics Challenge [56, 174-176] zeigt dennoch den 
enormen Fortschritt in der Kontaktmechanik der letzten Jahrzehnte. 


Im Rahmen von gleitenden Kontakten entsteht in der ersten Hälfte des 20. Jahrhunderts 
das Konzept der Blitztemperaturen durch Harmen Blok [31]. Darunter werden lokale, 
temporäre TemperaturerhGhungen verstanden, die deutlich über der Temperatur 
des Festkörpers liegen. Sie sind das Resultat von kurzzeitig in Berührung stehen- 
den Asperiten und treten bereits bei geringen Gleitgeschwindigkeiten auf. Neben 
Nachweisen in Experimenten [191, 223, 234] führen aufbauende Untersuchungen zur 
Temperaturentwicklung bei gleitenden Kontakten von Jaeger [114] und Archard [9] 
zu einer Serie von Veröffentlichungen, die sich weitestgehend mit der Abschätzung 
der mittleren oder maximalen Temperaturerhöhung bei bewegten oder stehenden 
Wärmequellen mit vordefinierter Ausdehnung auf der Oberfläche eines Halbraums 
beschäftigen [19, 65, 237]. Auch auf diesem Forschungsgebiet führt die verfügbare 
Rechenleistung Ende des 20. Jahrhunderts zu einer ganzen Reihe von Arbeiten, die 
unter anderem die Wärmeaufteilung [34, 245], den konvektiven Wärmeaustausch [149], 
die Wärmeentwicklung in beschichteten Körpern [205, 216, 226, 235], die experimen- 
telle Validierung [55] und die Auswirkungen räumlich verteilter Wärmequellen [4, 
54, 151, 243] untersuchen. Die Veröffentlichungen beziehen sich dabei zumeist auf 
die Sammlung von Wärmeleitungsproblemen von Carslaw und Jaeger [47]. Die 
Forschungsbemühungen zur Temperaturentwicklung werden in der Regel eher mit 
dem Verschleiß der Oberflächen als mit den Auswirkungen auf den Reibwert motiviert. 
Jedoch zeigen diverse Publikationen sowohl einen Einfluss der Temperatur auf die 
Materialparameter als auch auf die tatsächliche Kontaktfläche [17, 21, 44, 87, 157, 
181, 262] und damit sehr wohl einen Zusammenhang der Temperatur zur Reibkraft 
selbst [51, 171]. Der Einfluss auf die tatsächliche Kontaktfläche wird zum einen 
durch die direkte Beeinflussung der Materialparameter und zum anderen durch 
thermoelastische Verformungen der Oberflächenhöhen verursacht. Diese sind im 
Allgemeinen eng mit der Temperaturentwicklung verbunden. Das Phänomen der 
thermoelastischen Instabilität wird von James Richard Barber [16, 18, 20] sowohl 
experimentell als auch analytisch beschrieben und von Dow, Burton und Dundurs [39, 
40,67, 69] weiter untersucht. Nominell glatte Oberflächen treten hierbei an den höchsten 
Asperiten in Kontakt, die Reibungswärme führt zu einer Ausdehnung der ohnehin 
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schon in Kontakt stehenden Erhebungen. In der Folge erfahren diese eine ansteigende 
Belastung und damit einen wachsenden Wärmeeintrag. Dieser Kreislauf setzt sich fort 
und ist, mittlerweile unter der Bezeichnung thermoelastic dynamic instability (TEDI), 
immer noch Gegenstand aktueller Forschungen [1, 141]. Es findet zunehmend eine 
Weiterentwicklung von stationären Modellierungen zu transienten Lastfällen statt und 
grundlegende Kenntnisse über thermomechanische Einflüsse und deren Bedeutung 
für Reibkontakte entstehen [22, 32, 51, 109, 122, 171]. Neben den bereits genannten 
Auswirkungen wird auch der Einfluss der Temperatur und der Gleitgeschwindigkeit 
auf die Versagensscherspannung der Kontaktkörper erkannt. Veröffentlichungen aus 
den letzten Jahren beschäftigen sich unter Verwendung numerischer Hilfsmittel 
mit thermoelastischen Verformungen [60, 142, 146, 148] und Spannungen [147] von 
rauen Körpern, sowie mit Untersuchungen zu beschichteten Körpern [225, 264] und 
zu konvektivem Wärmeaustausch [163]. Außerdem entstehen Arbeiten zu thermo- 
mechanischen Kontakten [35, 49, 50, 140, 143, 224, 246], die die verschiedenen, 
bereits erwähnten mechanischen und thermischen Einflüsse in unterschiedlicher 
Ausprägung kombinieren. Grundsätzlich lassen sich dabei aufgrund der wesentlich 
komplexeren Zusammenhänge im Vergleich zu rein mechanischen Modellierungen 
deutlich weniger Veröffentlichungen finden, die jedoch ebenfalls eine konsequente 
Modellierung im Sinne einer transienten äquivalenten rauen Oberfläche vermissen 
lassen. Außerdem wird in dem Großteil der Arbeiten ebenfalls auf Mikroebene mit 
einem Coulomb’schen Reibgesetz gearbeitet. Inwieweit dies physikalisch sinnvoll ist, 
wird kritisch hinterfragt [75]. Die Komplexität der Modellierung und Simulation von 
tribologischen Kontakten wird darüber hinaus auch durch den Umfang aktueller 
Literaturübersichten [86, 162, 241, 250, 251] deutlich. 


Parallel zu den Forschungsbemühungen auf Mikroebene findet in der zweiten Hälfte 
des 20. Jahrhunderts eine Entwicklung von zwar phänomenologisch motivierten [2, 164, 
199, 200, 236], aber dennoch rein empirischen Reibgesetzen statt. Diese greifen häufig 
auf die Vorstellungen Coulombs von rauen Oberflächen als Ansammlung verformbarer 
Borsten zurück, weshalb sie des Öfteren unter der Bezeichnung Bürstenmodelle zusam- 
mengefasst werden. Grundsätzlich kann hierbei zwischen statischen und dynamischen 
Reibgesetzen unterschieden werden, wobei sich Letztere durch die Nutzung einer 
inneren Variable als Differentialgleichung schreiben lassen. Die auf diesem Wege 
sehr weit entwickelten Reibgesetze beinhalten zumeist sowohl eine Parametrisierung 
der Stribeck-Kurve im Sinne eines stationären, aber von der Relativgeschwindigkeit 
abhängigen Grenzfalls als auch eine Modellierung des partiellen Gleitens [197]. 

Aufgrund der Vielzahl von Reibgesetzen sollen nur einige ausgewählte dynamische 
Reibgesetze andiskutiert werden, um grundsätzliche Parallelen beider Forschungs- 
gebiete aufzuzeigen. Interessant ist beispielsweise der Vergleich des von Dahl [62] 
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erdachten dynamischen Reibmodells mit den analytischen Ergebnissen des partiell 
gleitenden Hertz-Kontaktes von Mindlin [169]. Dahl motiviert sein Modell durch die 
Analogie zu ideal elastisch-plastischem Materialverhalten. Die zugrundeliegenden 
Gleichungen beider Modelle lassen sich bei konstanter Last strukturell ineinander 
überführen [203]. Außerdem modelliert Dahl einen zusätzlichen stochastischen Anteil 
des Reibwerts als weißes Rauschen, ein Aspekt, der in weiterführenden dynami- 
schen Reibgesetzen nicht mehr aufgegriffen wird. Als Beispiele seien das LuGre- 
Modell [13], das elasto-plastische-Modell [71] und die Reibmodelle von Dietrich sowie 
Ruina [211] genannt. Die Vernachlässigung steht im Widerspruch zu Beobachtungen 
in Experimenten, die den stochastischen Charakter des Reibwerts immer wieder 
nahelegen [26, 66, 78, 133, 137, 228-230, 265]. Eine systematische Auswertung der 
experimentellen Daten hinsichtlich qualitativer Eigenschaften des stochastischen 
Reibwertanteils oder die gezielte Weiterentwicklung eines stochastischen Reibgesetzes 
wurden in der Literatur jedoch nicht gefunden. Frühe Charakterisierungsansätze von 
Rabinowicz [214, 215] stellen aber den Zusammenhang der Reibwertschwankungen 
zu den Oberflächenrauigkeiten her. Und auch neuere Messungen des Reibwerts 
beziehen die Oberflächenrauheit gezielt mit ein [158]. Erwähnenswert sind zudem 
auch Modellierungen mit temperaturabhängigen Reibwerten [159, 160, 182] sowie 
Veröffentlichungen zur Dynamik von Reibschwingern unter Berücksichtigung von 
thermoelastischen Ausdehnungen [212], die die tatsächlichen physikalischen Phäno- 
mene zwar stark abstrahieren, aber trotzdem den skalenübergreifenden Charakter von 
Systemen mit Reibung und die Bedeutung einer thermomechanischen Modellierung 
unterstreichen. Auf diesem Forschungsgebiet existieren ebenfalls zahlreiche und 
umfassende Zusammenfassungen in der Literatur [27, 77, 109, 110, 186, 187], die 
neben der Modellierung des Reibwerts auch auf ihre Folgen eingehen. 


Trotz fehlender Weiterentwicklung von stochastischen Reibgesetzen wird die Oberflä- 
chenrauigkeit immer wieder in der Modellbildung für dynamische Untersuchungen 
berücksichtigt. Neben Publikationen zu reinen Vertikalschwingungen von rauen 
Körpern [93-97, 115, 230, 258], zu Schwingungen von Mehrfreiheitsgradsystemen mit 
rauen Oberflächen und deterministischen Reibwerten [78, 220, 228] sowie berechneten 
Reibwerten [25] lassen sich auch Veröffentlichungen zu rein tangentialen Freiheitsgra- 
den mit stochastischem Reibwert im Sinne eines unphysikalischen weißen Rauschens 
finden [79, 83, 100]. Forschungsbemühungen zu Mehrfreiheitsgradsystemen unter 
Berücksichtigung eines stochastischen Reibwerts im Sinne eines farbigen und/oder 
korrelierten Rauschprozess wurden in diesem Zusammenhang nicht gefunden. 


Zusammenfassend kann festgehalten werden, dass das interdisziplinäre und mehrere 
Größenskalen überspannende Forschungsfeld der Tribologie in vielen Bereichen bereits 
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gut erschlossen ist. Die tiberwiegende Mehrheit der Arbeiten konzentriert sich jedoch 
auf den quasi-statischen Normalkontakt. Neben der Vielzahl an Einflussfaktoren treten 
bei gleitenden Kontakten im Vergleich zu statischen Kontakten zwei Besonderheiten 
auf. Erstens kann für eine konsequente Modellierung zweier gleitender Körper mit 
rauen Oberflächen nicht auf eine statische, äquivalente raue Oberfläche zurückgegriffen 
werden. Der aktuelle Stand der Forschung weist diesbezüglich noch eine Lücke auf. 
Zweitens werden die Wärmeentwicklung und ihre Folgen für den gleitenden Kontakt 
von Bedeutung sein. Die Auswirkungen der Temperaturentwicklung auf die Materi- 
alparameter und den thermomechanischen Kontakt sind zumindest in der Theorie 
und aus Experimenten bekannt, detaillierte Simulationsmodelle fehlen jedoch. Um 
die Einflüsse beider Aspekte auf den Reibwert abzubilden, ist die Verwendung eines 
vorgegebenen Reibwerts auf Mikroebene zu vermeiden. 

Des Weiteren fehlen über das weiße Rauschen hinausgehende Modellierungen des 
stochastischen Reibwertanteils. Dies liegt neben der nicht vorhandenen Weiterent- 
wicklung von stochastischen Reibgesetzen auch an nicht zielführenden Auswertungen 
vorhandener Experimente und fehlenden Kontaktsimulationen mit einer konsequenten 
Modellierung der zweiseitigen Rauigkeit. Eine entsprechende Lücke des Forschungs- 
standes lässt sich dementsprechend bei weiterführenden Studien zu den Auswirkungen 
von stochastischen Reibwerten auf die Dynamik von technischen Systemen verorten. 


1.3 Struktur und Zielsetzung der Arbeit 


Nach der erfolgten Motivation und der Zusammenfassung des aktuellen Forschungs- 
standes schließt eine Beschreibung der Struktur und der Zielsetzung der Arbeit das 
Kapitel 1 ab. Dafür werden zunächst die Ziele der Arbeit im Sinne des folgenden, 
offenen Forschungsbedarfes abgeleitet: 


e Die Berechnung des Reibwerts für zwei gleitende Körper mit rauen Oberflächen 
bei unterschiedlichen Relativgeschwindigkeiten und Lasten. 


e Die konsequente Modellierung der Oberflächenrauigkeiten von beiden Kontakt- 
körpern. 


Die Berücksichtigung von Reibungswärme und thermoelastischen Verformungen. 


Die physikalisch sinnvolle Modellierung und Parametrisierung des stochastischen 
Reibwertanteils. 


e Die Untersuchung der Auswirkungen etwaiger stochastischer Anteile des Reib- 
werts im Kontext von reibungserregten Schwingungsphänomenen. 


10 


1.3 Struktur und Zielsetzung der Arbeit 


Wahrend diese Forschungslticken geschlossen werden, soll zudem der Aspekt einer 
anschließenden Validierungsmöglichkeit der Simulationsergebnisse durch Experimente 
sichergestellt sein. Die Oberflächenhöhen im Kontaktbereich können vermessen werden, 
sodass eine diskrete Menge an x, y und z-Werten zur Verfügung steht. Die mittlerweile 
von modernen Messgeräten erreichbaren Auflösungen von 128 x 128 bis 2048 x 2048 
Datenpunkten müssen deshalb als Eingangsgrößen für das Simulationsmodell ver- 
wendet werden können. Das Durchführen der Simulationen für Rechengitter in dieser 
Größenordnung innerhalb vertretbarer Rechenzeiten, auch auf handelsüblichen Com- 
putern, kann dabei als untergeordnetes Forschungsziel formuliert werden. 

Die Arbeit ist wie folgt strukturiert. Kapitel 2 legt zunächst die Grundlagen für die 
weitere Arbeit. Dazu gehören generelle Aspekte zur Modellierung und Simulation von 
rauen Oberflächen und Reibwerten in Unterkapitel 2.1. Verschiedene Charakterisie- 
rungsmöglichkeiten von stochastischen Prozessen werden in Abschnitt 2.1.1 aufgezeigt 
und für raue Oberflächen in Abschnitt 2.1.2 und für Reibwerte in Abschnitt 2.1.3 
diskutiert. Die beiden letztgenannten Abschnitte beinhalten darüber hinaus mögliche 
Simulationsverfahren für die Erstellung von rauen Oberflächen und von Reibwerten 
mit gewünschten Eigenschaften. In Unterkapitel 2.2 schließt sich die Beschreibung von 
Kontinua im Rahmen eines linear thermoelastischen Materialverhaltens an, die sich in 
eine Beschreibung unterschiedlicher Lösungsmethoden für die erhaltenen, partiellen 
Differentialgleichungen in Abschnitt 2.2.1 sowie ausgewählter numerischer Methoden 
in Abschnitt 2.2.2 gliedert. Ein Überblick zu reibungserregten Schwingungen wird in 
Unterkapitel 2.3 gegeben. Dafür werden die zugrundeliegenden Mechanismen anhand 
zweier Minimalmodelle erörtert. Hierzu zählt einerseits der negative Reibwertgradient 
in Abschnitt 2.3.1 und andererseits die nicht-konservative Kopplung in Abschnitt 2.3.2. 
Beide Minimalmodelle repräsentieren eine Vielzahl von technischen Systemen. 

Die gelegten Grundlagen dienen anschließend in Kapitel 3 für die thermomechanische 
Modellierung zweier gleitender Körper mitrauen Oberflächen. Dazu zählt zunächst eine 
grundsätzliche Beschreibung der getroffenen Annahmen, bevor die zugehörigen parti- 
ellen Differentialgleichungen in Unterkapitel 3.1 für mechanische Oberflächenlasten in 
Abschnitt 3.1.1, für thermische Oberflächenlasten in Abschnitt 3.1.2 und schließlich für 
allgemeine Oberflächenlasten in Abschnitt 3.1.3 gelöst werden. Unterkapitel 3.2 beinhal- 
tet neben einer Verifikation der hergeleiteten Lösungen für verschiedene mechanische 
und thermische Lastfälle durch analytische Vergleichslösungen aus der Literatur auch 
eine Diskussion der erhaltenen Ergebnisse und ausgewählter, weiterführender Aspekte. 
Unterkapitel 3.3 vervollständigt schlussendlich die Modellbildung und beschreibt den 
Simulationsablauf. 

Die Ergebnisse der thermomechanischen Kontaktsimulationen werden in Kapitel 4 
für verschiedene Anwendungsfälle vorgestellt und im Sinne des offenen Forschungs- 
bedarfs ausgewertet. Dazu gehören einerseits Kontaktkonfigurationen mit isotropen 
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Oberflachenstrukturen in Unterkapitel 4.1 und andererseits Kontaktkonfigurationen mit 
anisotropen Oberflächenstrukturen in Unterkapitel 4.2. Es schließt sich eine Zusammen- 
fassung und ein Vergleich der Ergebnisse in Unterkapitel 4.3 an. Des Weiteren werden 
die Zusammenhänge zwischen charakteristischer Oberflacheneigenschaften und den 
erhaltenen kontaktmechanischen Größen herausgearbeitet, bevor die Erweiterung um 
temperaturabhängige Materialparameter untersucht wird. 

Die erhaltenen Ergebnisse werden anschließend in Kapitel 5 in den Kontext von 
reibungserregten Schwingungen gestellt. Dies gelingt durch eine Berücksichtigung 
des berechneten Reibwerts in den Minimalmodellen aus Unterkapitel 2.3, wobei sich 
die Diskussion in Unterkapitel 5.1 auf den deterministischen und in Unterkapitel 5.2 
auf den stochastischen Reibwertanteil konzentriert. Für Letzteren erfolgt zunächst 
die im Sinne des offenen Forschungsbedarfs liegende Modellierung des stochastischen 
Reibwertanteils in Abschnitt 5.2.1, bevor die Auswirkungen auf reibungserregte Schwin- 
gungen in Abschnitt 5.2.2 herausgearbeitet werden. Beide Unterkapitel beinhalten eine 
Stabilitätsanalyse des stationären Betriebspunktes durch jeweils geeignete Methoden. 
Die Ergebnisse werden insbesondere vor dem Hintergrund bereits bekannter oder 
eventuell neuer Instabilitätsmechanismen bei technischen Systemen mit trockener 
Reibung ausgewertet. 

Die Arbeit schließt in Kapitel 6 mit einer Zusammenfassung und einem Ausblick ab. 
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Dieses Kapitel setzt die Grundlagen fiir die weiteren Ausfiihrungen. Es kniipft dabei 
an den Stand der Forschung an und zeigt zunächst in Unterkapitel 2.1 verschiedene 
Modellierungs- und Simulationsmöglichkeiten für raue Oberflächen und Reibwerte 
auf. Die der Kontaktmechanik zugrundeliegenden Gleichungen und ausgewählte 
numerische Aspekte werden anschließend in Unterkapitel 2.2 diskutiert, bevor ein 
Überblick über reibungserregte Schwingungen in Unterkapitel 2.3 gegeben wird. 


2.1 Modellierung und Simulation von rauen 
Oberflächen und Reibwerten 


Sowohl raue Oberflächen als auch Reibwerte können als stochastischer Prozess Xs 
aufgefasst werden. Unter einem stochastischen Prozess versteht sich dabei eine in 
s stetige Folge zufälliger Vorgänge. Die unabhängigen Parameter s = [x,y,...,t]" 
repräsentieren beispielsweise den Ort und/oder die Zeit. Da sowohl raue Oberflächen 
als auch Reibwerte zentraler Bestandteil dieser Arbeit sind, werden die wichtigsten Cha- 
rakterisierungsmöglichkeiten von stochastischen Prozessen zunächst in Abschnitt 2.1.1 
genauer erläutert. Anschließend findet eine entsprechende Einordnung von rauen 
Oberflächen in Abschnitt 2.1.2 sowie von Reibwerten in Abschnitt 2.1.3 statt, wobei 
mögliche Simulationsmethoden beschrieben werden. 


2.1.1 Stochastische Prozesse 


Ausgangspunkt sei ein stochastischer Prozess X, mit unabhängigem Parameter s und 
zugehöriger Verteilungsdichte px(x,s). Die folgenden Ausführungen orientieren sich 
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an Natke [177] und lassen sich auf beliebig viele unabhangige Parameter s tibertragen. 
Die Verteilungsdichte px(x,s) lässt sich durch die Momente i-ter Ordnung 


+00 
mË (s) SE [xi] = J x'px(x, s)dx 


(0) 


charakterisieren, das Moment nullter Ordnung m X 


= 1 legt die Normierung fest. Die 
zentralen Momente i-ter Ordnung sind mit 


u® (s) =E (x — mis) | = fk - mo) px(x,s)dx 


gegeben. Eine Aussage über den Zusammenhang zweier Zufallsvariablen X;, und Xs, 
kann die Autokovarianzfunktion 


Kxx(s1,82) =E | (Xs, - m\(s1)) (Xa well 


+00 +00 


St d (x = my? (s1)) (v = EN px(x, y, $1, $2)dxdy 


—00 —00 


mit der Verbundwahrscheinlichkeitsdichte px(x, y, $1, 52) liefern. 

Im Folgenden soll eine Beschränkung auf stationäre, mittelwertfreie Prozesse erfolgen, 
für die px(x,s) = px(x) und mË 
anzfunktion Kxx(s1, 52) in die Autokorrelationsfunktion Rxx(As) über und lässt sich 


= 0 gilt. Für solche Prozesse geht die Autokovari- 


durch eine Normierung auf das zentrale Moment zweiter Ordnung u im Intervall 
[-1,1] darstellen. Sie ist außerdem eine gerade Funktion und hängt nur noch von der 
Differenz As = |sı - s2| ab. 

Der Zusammenhang zwischen Autokorrelationsfunktion Rxx(As) und zweiseitigem 
Leistungsdichtespektrum Sxx(w) ist durch das Wiener-Chintschin-Theorem und die 
kontinuierliche Fouriertransformation 


F(f(s)}w) = f f(s) Tds (2.1) 
FUN) =z | fiel" 22 


oo 
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als 
Sxx(w) = F {Rxx}(o) 


gegeben. Der Zusammenhang zwischen zweiseitigem Leistungsdichtespektrum und 
dem zentralen Moment zweiter Ordnung lässt sich über 


+00 


1 
Rxx(As = 0) = al Sxx(w)dw = ike (2.3) 


00 


herstellen. Da das zweiseitige Leistungsdichtespektrum eine gerade Funktion und für 
reelle Signale positiv ist, wird es auch oft als einseitiges Spektrum zusammengefasst. 
Bei ergodischen Prozessen ist die Kenntnis der Verteilungsdichte px(x) nicht mehr 
Voraussetzung für die Berechnung der Erwartungswerte, es reicht die Kenntnis einer 
einzelnen Realisierung x) (s) aus 


+Ts 
mÊ = Jo = (x (s)] ds 
-T 
+Ts 
HS = tim. ae | (emp!) as, 
-T, 
AT. 


1 ; , 
Kxx(As) = lim FA A (x (s) - my’) (x (s + As) - my’) ds. 
gre. S 
E 


Da in der Praxis nur eine begrenzte Anzahl an endlichen Realisierungen x) (s) = 
x(mAs) mit m = 0,...,M -1,As = T,/M und j = 1,...,J zur Verfügung steht, 
müssen die eingeführten Kenngrößen geschätzt werden. Für die in dieser Arbeit 
relevanten Momente ergeben sich die erwartungstreuen und konsistenten Schätzungen 
für Mittelwert 


1 M-1 
i” — Ds x)(mAs) 
m=0 
Varianz 
M-1 
1 2 
TA = MI A (xmas) — mQ) z 
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Schiefe 
M-1 3 
x (mAs) — m) 
Ex = M-IM-2 3 
LO 
Hx 


und Kurtosis 


E (xmas) - i) 

~(4) _ M(M + 1) m=0 S (M = De +3 

Ex = (Mf —1)(M — 2)(M -3) Di (M-2)M-3) ~ 
x 


Für die Autokorrelationsfunktion ergibt sich bei Voraussetzung eines periodischen 
Signals mit x(mAs) = x)((m + M)As) die für ergodische Prozesse erwartungstreue, 
allerdings nur für normalverteilte Prozesse konsistente Schätzung 


M-1 
mae) 1 ; ; 
RY (KAS) = se A, xM(mAs)x (m+ k)As), ` Ech. Mi 


m=0 


Das zweiseitige Leistungsdichtespektrum lässt sich zunächst als Produkt der komplex- 
konjugierten, finiten Fouriertransformierten einer Realisierung des Prozesses x) (s) 
T, 2 


1 ’ 1 
Sxx(w) = lim z E [IF {1x0 Ho, T] = lim zE yA xM(s)eT®sds 
sco E Dk ASA 


ausdriicken, woraus mit der diskreten Fouriertransformation 


M-1 
DF {f(mAs)} Mag (or) =As A. f(mAs)e "M,  k=0,...,M-1 (24) 


m=0 


M-1 
DF (ta) (mAs) =—— Y Plate ann,  m=0,...,M-1 (25) 
MAs el 


und wx = Zmk/M die asymptotisch erwartungstreue, aber nicht konsistente Schätzung 


1 


ans PF EP Hor) / 


Sxx(j)(@x) = 
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folgt. Sie wird auch Periodogramm genannt und steht über die diskrete Fouriertrans- 
formation in direktem Zusammenhang zur Autokorrelationsfunktion 


ED (wx) = DF RY Mow), k=0,...,M-1. 


Da hierbei die Erwartungswertbildung entfällt, kann die Schätzung durch anschlie- 
ßende Mittelung über mehrere Realisierungen häufig noch verbessert werden. 


2.1.2 Raue Oberflächen 


Technische Oberflächen sind grundsätzlich gewissen Fertigungstoleranzen und damit 
Gestaltabweichungen unterworfen. Der Begriff Rauigkeit kann dabei durch die Norm 
DIN 4760 als Gestaltabweichung dritter bis fünfter Ordnung identifiziert werden und 
lässt sich damit zwischen den Gestaltabweichungen Welligkeit und Gitteraufbau des 
Werkstoffs einordnen. Der Übergang zwischen den einzelnen Größenordnungen ist 
hierbei fließend und kann je nach Anwendungsfall unterschiedlich sein. Für die Vermes- 
sung von Oberflächen stehen mittlerweile zahlreiche berührende und berührungslose 
Messmethoden zur Verfügung [244]. Diese liefern eine diskrete Menge an x, y, und 
z-Werten, wobei für eine möglichst umfassende Beschreibung der Gestaltabweichungen 
in der Regel mehrere Messungen mit unterschiedlichen Auflösungen kombiniert werden 
müssen [111, 202]. 

Werden die erhaltenen Höhenwerte z als Realisierung eines ortsabhängigen, statio- 
nären und ergodischen Prozesses Zs, s = [x,y]' mit der Verteilungsdichtefunk- 
tion pz(z) aufgefasst, können für die Beschreibung des Höhenprofils die bereits in 


Abschnitt 2.1.1 aufgeführten Methoden herangezogen werden. In diesem Kontext 
(1) 
Z y 


Varianz R3 = uA als Quadrat der Standardabweichung R4, Schiefe Rsk = a IR: 
und Kurtosis Rx, = Msg / R3 durch die Norm DIN EN ISO 4287 etabliert, während die 
entsprechenden Bezeichnungen Sm, Sq, Ssk und Sku im Rahmen von Flächenmessungen 
durch die Norm DIN EN ISO 25178 verwendet werden. Die Literatur beinhaltet 
viele weitere Parameter, deren Aussagekraft jedoch in den meisten Fällen zumindest 


haben sich im Rahmen von Linienmessungen die Bezeichnungen Mittelwert Rm = m 


im Zusammenhang mit kontaktmechanischen Fragestellungen umstritten ist. Eine 
entsprechende Sensitivitätsanalyse wird unter anderen von Duo et al. [70] durchgeführt. 
Vermessungen von neuen Oberflächen weisen zudem häufig normalverteilte Oberflä- 
chenhöhen Zs ~ N (Sm, si auf [88, 178], weshalb viele Forschungsbemtihungen auf 
einer solchen Verteilung basieren. Jedoch unterliegen die Oberflachen, beispielsweise 
aufgrund von Verschleiß, ständigen Veränderungen. Die Modellierung mit normal- 
verteilten Oberflächenhöhen erscheint daher für die praktische Anwendung selten 
gerechtfertigt. Darüber hinaus ist die Annahme eines stationären und ergodischen Pro- 
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zesses, wenn überhaupt, nur für stationäre Betriebsbedingungen und den eingelaufenen 
Zustand möglich. Eine Auflistung von typischen Kennwerten technischer Oberflächen, 
die aufgrund verschiedener Fertigungsverfahren resultieren, lässt sich bei Gomeringer 
et al. [84] finden. 

Neben den charakteristischen Kenngrößen der Verteilungsdichte liefern Autokorrela- 
tionsfunktion Rxx(Ax, Ay) oder Leistungsdichtespektrum $zz (wx, wy) weitere Infor- 
mationen über die Eigenschaften der Oberflächen. Dabei hat sich mittlerweile die 
in Abbildung 2.1 gezeigte Parametrisierung des Leistungsdichtespektrums für den 
isotropen Fall in der Form 


0 für læll < a, 
a H für on < ||w|| < wr, 

Szz(w) = S04 an 2.6) 
lo für w, < lol < ws, 


0 für ws < ||@|| 


etabliert [111, 202]. Diese beinhaltet neben den Kreisfrequenzen ||w|| = Jo? + Wy und 
dem Hurst-Exponenten H = 3 — Dy als Funktion der fraktalen Dimension Dy auch 
die Konstante So, die eine beliebige Normierung der Standardabweichung S4 gemäß 
Gleichung (2.3) zulässt. Eine entsprechende Verallgemeinerung auf den anisotropen 


Fall mit ||@|| = zeg + drng, dx # Py oder beliebigen weiteren Parametrisierungen 
des Leistungsdichtespektrums ist leicht möglich. Bei Vorgabe eines Gitters sind die 
in Gleichung (2.6) enthaltenen Kreisfrequenzen bereits durch die Diskretisierung 
vorgegeben. Des Weiteren lassen Vermessungen von technischen Oberflächen die 
Eingrenzung des Hurst-Exponenten auf den Bereich H ~ 0.8 + 0.1 zu [201, 202]. 


ka 

= 

P 

oN 
BN 

n -2(1+H) 

Ke 

S 

w] Wy Ws 


log ||w|| in m! 


Abbildung 2.1: Leistungsdichtespektrum > (w) gemäß Gleichung (2.6) 


Für kontaktmechanische Simulationen ist es wegen der zeit- und materialaufwendigen 
sowie skalenübergreifenden Vermessungen von Oberflächen sinnvoll, auf künstlich 
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(a) Szz (w) aus Gleichung (2.6) mit (b) Szz (w) aus Gleichung (2.6) mit (c) Szz (w) aus Gleichung (2.6) mit 
wr = 80m mm” wr, =8n mm! wr = 0.87 mm” 


(d) Szz (w) aus Gleichung (2.6) mit (e) Szz (w) aus Gleichung (2.6) mit (f) Szz (w) aus Gleichung (2.6) mit 
wr = 8r mm! und Ssk = —0.7, wr = 8n mmh, d = 1, Yy =3 wr = Bot mm, px = 5, Yy =1 
Sku =3 und Sk = —0.7, Sku =3 

Abbildung 2.2: Verschiedene Realisierungen z0) (x,y), j =1,...,6 rauer Oberflächen der Länge Ly = Ly = 

1 mm mit Hurst-Exponent H = 0.8 - wenn nicht anderweitig angegeben gilt Zs ~ N (0, 2) 


erzeugte Oberflächen mit möglichst gleichen Eigenschaften zurückzugreifen. Für 
diesen Zweck existiert mittlerweile eine ganze Reihe verschiedener Verfahren für 
die Erzeugung rauer Oberflächen [3, 194]. Welche Methode letztendlich verwendet 
wird, hängt von den spezifischen Anforderungen ab. Das für diese Arbeit verwendete 
Simulationsprogramm baut auf den Veröffentlichungen von Schreiber und Schmitz [222] 
sowie Perez und Almquvist [198] auf. Die Vorgehensweise erlaubt die Erzeugung 
rauer Oberflächen bei gleichzeitiger Vorgabe des Leistungsdichtespektrums Szz (w) 
und der Verteilungsdichte pz(z). Beide Größen können auch aus bereits durchge- 
führten Vermessungen von Oberflächen stammen und lassen damit eine direkte 
Verknüpfung mit Experimenten zu. Verschiedene Realisierungen z0 (x, y), j = 1,...,6 
von rauen Oberflächen, die mit dieser Simulationsvorschrift für unterschiedliche 
Leistungsdichtespektren Szz (w) und Verteilungsdichten pz(z) erstellt wurden, sind in 
den Abbildungen 2.2(a)-(f) gezeigt. 

Für den in Abbildung 2.3 skizzierten Simulationsablauf wird eine bereits mit der 
geforderten Verteilungsdichte pz(z) vorliegende Menge an Höhenwerten zp, der Dimen- 
sion M x N durch das Programm iterativ so lange umsortiert, bis ebenfalls das 
gewünschte Leistungsdichtespektrum Szz(w) erreicht wird. Dies gelingt aufgrund 
der noch offenen Freiheitsgrade der Phasenverschiebungen. In der Regel wird das 
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Leistungsdichtespektrum 
Szz(wx, w y ) 


Solange Tell > eu berechne: 


Ies DF {zy,}wx, Wy) 


Realisierung z(x,y) = Zp, 
mit Verteilungsdichte pz(z) 


. Ppz = arctan (Im (c) /Re (c)) 


» 2877 = DF { VSzz exp LU oe (x,y) 


. Ersetze den i-ten Wert 


i 
< TEE < TE 
2577 Zëss e 


durch den i-ten Wert 


1 
Zpz re Zp <:i 
und erhalte so umsortiertes zp; (x, y) 


Abbildung 2.3: Simulationsablauf für die Erzeugung stationärer, ergodischer Oberflächen oder Reibwerte 


geforderte Leistungsdichtespektrum Szz(w) wegen der begrenzten Anzahl an Wer- 
ten nicht exakt erreicht, weshalb die Iterationen bei Erreichen einer vordefinierten 
Toleranz eut abgebrochen werden. Die so umsortierten Höhenwerte zp, besitzen die 
geforderte Verteilungsdichte pz(z) und erfüllen näherungsweise das Leistungsdich- 
tespektrum Szz(w), während die so erhaltenen Höhenwerte zs,, dem vorgegebenen 
Leistungsdichtespektrum Szz(w) folgen. 

Da das Erzeugen der Oberflächen kein aktiver Teil der Kontaktsimulation selbst ist, lässt 
sich die benötigte Rechenzeit dem Pre-Processing zuordnen und fällt nicht weiter ins 
Gewicht. Auf handelsüblichen Rechnern ermöglicht die Verwendung der schnellen Fou- 
riertransformation zudem eine Erzeugung der Oberflächen in anwendungsrelevanten 
Gittergrößen von bis zu 2048 x 2048 Elementen innerhalb weniger Minuten. 


2.1.3 Reibwerte 


Viele Messungen in der Literatur deuten auf den stochastischen Charakter von trockenen 
Reibwerten hin [26, 66, 78, 133, 137, 228-230, 265], sodass deren Modellierung als 
stochastischer Prozess M; mit zugehöriger Verteilungsdichte p (u, s) gerechtfertigt 
erscheint. Ganz analog zu den Vermessungen von Oberflächen wird auch hier nur eine 
beschränkte Anzahl an endlichen Realisierungen ut (s ) j =1,...,J erhalten, weshalb 
für eine Charakterisierung der erhaltenen Reibwertverläufe mithilfe der Methoden aus 
Abschnitt 2.1.1 zwangsläufig auf eine Vermessung während stationärer Betriebspunkte 
zu achten ist. Dadurch lässt sich auch eine generelle Einschränkung auf stationäre, 
ergodische Prozesse M,; mit zugehöriger Verteilungsdichte pm(u) motivieren. Für 
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dynamische Simulationen bietet sich dementsprechend zunächst die Darstellung des 
Reibwerts Ms = ua + us als Summe eines deterministischen Anteils ua = E[M;] und 
eines mittelwertfreien stochastischen Anteils us = Ms - E [Ms] an. 

Der deterministische Anteil ua lässt sich beispielsweise mit einer der in der Literatur 
zahlreich vorhandenen Parametrisierungen der Stribeck-Kurve für trockene Kontakte 
identifizieren [197]. Für eine explizite Parameteridentifikation sind in der Praxis zeit- 
und materialaufwendige Messungen nötig, da die Stribeck-Kurve einen stationären, 
aber von der Relativgeschwindigkeit abhängigen Grenzfall darstellt. Der stationäre Wert 
des Reibwerts wird dabei erst nach einer gewissen Zeitspanne beziehungsweise einem 
bestimmten Gleitweg erreicht. Abbildung 2.4 zeigt verschiedene qualitative Verläufe 
der Stribeck-Kurve als Funktion der relativen Gleitgeschwindigkeit. Viele Forschungs- 
bemühungen, insbesondere im Kontext von reibungserregten Schwingungen, beruhen 
auf geschwindigkeitsabhängigen Verläufen der Reibkraft, wobei zumeist der Gradient 
der Kurve von entscheidender Bedeutung ist. Eine Berechnung des Kurvenverlaufs 
ohne Messungen ist daher nicht nur aus Kostengründen, sondern auch im Sinne einer 
gezielten konstruktiven Beeinflussung des Reibwerts erstrebenswert. Aufgrund der 
vorhandenen Parametrisierungen bereiten analytische und simulative Behandlungen 
eines solchen funktionellen Zusammenhangs jedoch keine Probleme. 


Hay 


0 || rel in m/s 


Abbildung 2.4: Drei verschiedene qualitative Verläufe des deterministischen Reibwertanteils ua 


Anders sieht es mit einer simulativen Umsetzung des stochastischen Anteils u; aus. Die- 
ser kann zwar prinzipiell analog zu der in Abbildung 2.3 gezeigten und in Abschnitt 2.1.2 
beschriebenen Vorgehensweise für die Erzeugung von Oberflächenhöhen erstellt wer- 
den, hierbei fallen jedoch zwei Probleme auf. Zum einen liefert die Literatur im 
Gegensatz zu rauen Oberflächen wenig bis gar keine Informationen über die Ver- 
teilungsdichte pyu(u) oder das Leistungsdichtespektrum Smm(w) von Reibwerten. 
Eine Quantifizierung dieser beiden Größen - sei es durch die Berechnung eines 
stochastischen Reibwerts oder die Auswertung von Messungen - ist daher notwendig. 
Zum anderen fällt die Rechenzeit für die Erzeugung des Reibwerts aufgrund der 
wesentlich häufigeren Aufrufe bei Parameterstudien deutlich mehr ins Gewicht. Es 
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erscheint vorteilhafter, den Reibwert mit möglichst geringem Mehraufwand während 
der eigentlichen Simulation nebenher zu erzeugen. Diese Vorgehensweise ist bei 
deterministischen Reibwerten bereits weit verbreitet (vgl. dynamische Reibgesetze in 
Unterkapitel 1.2) [197]. Die bereits von Dahl [62] vorgeschlagene Modellierung us = Zs 
als normalverteiltes, weißes Rauschen ist jedoch aufgrund der unbeschränkten Varianz 
des weißen Rauschens unphysikalisch. Als physikalisch sinnvollere Modellierung 
bietet sich in diesem Zusammenhang aber ein farbiges Rauschen an. Entsprechende 
Simulationsansätze lassen sich vielfach in Verbindung mit Fahrbahnanregungen in der 
Fahrzeugdynamik finden [129, 209]. Hierbei wird das gewünschte Leistungsdichte- 
spektrum oder die gewünschte Autokorrelationsfunktion zunächst approximiert und 
anschließend durch Filtergleichungen umgesetzt. Diese lassen sich im mathematischen 
Sinne als System linearer, stochastischer Differentialgleichungen erster Ordnung 


X, = AX, + BS, Xa = Xo (2.7) 


mit Driftterm AX, und Diffusionsterm B darstellen, wobei die Anregung ©, ein nor- 
malverteiltes, weißes Rauschen darstellt. Aufgrund der Linearität von Gleichung (2.7) 
ist deren stationäre Lösung ebenfalls normalverteilt. Für nicht-normalverteilte Prozesse 
kann im Allgemeinen auf nichtlineare Drift- und Diffusionsterme 


X, = A(X,) + B(X,) Bs, Xs = Xo (2.8) 


zurückgegriffen werden. Entsprechende Beispiele lassen sich in [57, 130, 227, 266] 
finden. Grundsätzlich erfolgt für die Simulation aber in allen Fällen zunächst die 
formale Überführung in eine Interpretation nach Itö 

1 OB(X;s) 


dX, = [a (Xs) +5 ze 


B(X.)| ds + B(X;)dW,, Xs, = Xo, (2.9) 


wobei das Wiener Inkrement dW, = 2ds den Zusammenhang zum weißen Rau- 
schen £; abbildet. Die Interpretation nach Itö hat unter anderem den Vorteil, dass 


s2 
E [Bean =0 


S1 


gilt. Es sei ergänzt, dass sich für das stochastische Integral in Gleichung (2.8) bei 
einer anderen Wahl der Stützstellen ebenso die Interpretation nach Stratonovich 
finden lässt. Durch entsprechende Korrekturterme können die beiden Interpretationen 
jedoch immer ineinander überführen werden [11]. Da die Lösung der stochastischen 
Differentialgleichung ebenfalls ein stochastischer Prozess ist, lässt sich dieser wiederum 
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mit den Methoden aus Abschnitt 2.1.1 charakterisieren. Wird der stochastische Anteil 
des Reibwerts beispielsweise durch us = Xıs identifiziert, lässt sich die zugehörige 
Randverteilungsdichte pu(u) durch 


+00 +00 
pmd = fo f pxladxa...dx 


aus der Verbundwahrscheinlichkeitsdichtefunktion px(x) berechnen. Drift- und Dif- 
fusionsterm in Gleichung (2.9) werden im späteren Verlauf der Arbeit so gewählt, 
dass sowohl die gewünschte Verteilungsdichte p m(u) als auch das gewünschte Leis- 
tungsdichtespektrum S um (w) erhalten werden. Die dafür notwendigen Informationen 
liefern thermomechanische Kontaktsimulationen. An dieser Stelle ist noch einmal 
hervorgehoben, dass die stochastische Differentialgleichung (2.8) dafür eine stationäre 
Lösung besitzen muss. 

Für das Lösen der Gleichung (2.9) bieten sich neben den numerischen Integrations- 
verfahren für stochastische Differentialgleichungen auch die Evolutionsgleichungen 
für die Verbundwahrscheinlichkeitsdichtefunktion in der Form von partiellen Diffe- 
rentialgleichungen an. Für eine ausführlichere Diskussion von numerischen Aspekten 
wird auf die Literatur verwiesen [14, 128, 134]. Da für dynamische Simulationen aber 
in der Regel immer nur einzelne Realisierungen u) (s) von Interesse sind, kann das 
eigentliche Differentialgleichungssystem des physikalischen Systems um das Glei- 
chungssystem (2.9) erweitert und die Realisierung u) (s) während der eigentlichen 
Simulation erzeugt werden. Dafür wurde ein Runge-Kutta-Verfahren für stochastische 
Differentialgleichungen implementiert [218, 219]. Dieses hat unter anderem den Vorteil, 
dass es im Vergleich zum oft verwendeten Euler-Maruyama-Verfahren auch eine höhere 
Konvergenzordnung bezüglich des deterministischen Integrals besitzt. 


2.2 Thermoelastischer Festkörper 


Für die mathematische Beschreibung von Kontinua stehen fünf universell geltende 
Bilanzgleichungen (Masse-, Impuls- und Drehimpulsbilanz sowie erster und zweiter 
Hauptsatz der Thermodynamik) zur Verfügung, aus denen sich mithilfe von Materi- 
algleichungen partielle Differentialgleichungen für die Verschiebungs- und Tempera- 
turfelder herleiten lassen. In Abschnitt 2.2.1 werden die so erhaltenen Gleichungen für 
hyperelastische Materialien und lineare Thermoelastizität erläutert, bevor ausgewählte 
numerische Aspekte in Abschnitt 2.2.2 folgen. 
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2.2.1 Lineare Thermoelastizitat 


Ein homogener, isotroper Festkörper kann im Rahmen der linearen Thermoelastizität 
durch die beiden linearen, gekoppelten, partiellen Differentialgleichungen für das 
Verschiebungsfeld r = (u vw)" mit Vr < 1 und für das Temperaturfeld 0 = T - Tù als 
Differenz der absoluten Temperatur T zu einer Referenztemperatur To mit |O/To| « 1 


uAr + (A +u)V(V-r) + pf =A +2u) aVO + pDyo, (2.10) 


ð 98 
-V -q - (8A +2) aV -r =pcas -Q (2.11) 


beschrieben werden [33]. Für deren grundlegende Herleitung aus den thermomechani- 
schen Bilanz- und Materialgleichungen wird auf entsprechende Literatur verwiesen [33, 
184, 185, 193, 239], an der sich auch die folgenden Ausführungen orientieren. Die 
Bilanzgleichungen werden durch die Materialgleichungen von Hooke, Fourier und 
Gibbs 


o =u (Vr +(Vr)') + (AV -r — (3A + 2u) a0) I, (2.12) 
q =- KVO, (2.13) 
-fig + CAPA oyy (2.14) 

0 


ergänzt, die den Zusammenhang des Temperatur- und Verschiebungsfeldes zum Span- 
nungstensor o, den Zusammenhang des Temperaturfeldes zum Wärmestrom q und den 
Zusammenhang des Temperatur- und Verschiebungsfeldes zur spezifischen Entropie s 
festlegen. Die Gleichungen (2.10)-(2.14) beinhalten neben der zeitlichen Änderung 
des räumlichen Geschwindigkeitsfeldes D;v, der äußeren Volumenkraftdichte f und 
Wärmestromdichte Q auch die beiden Lamé-Konstanten u = E/2(1+v) und A = 
Ev/(1+v)(1-2v) als Funktionen des Elastizitätsmoduls E und der Querkontraktions- 
zahl v, den Wärmeausdehnungskoeffizient a sowie den Wärmeleitfähigkeitstensor K, 
die Dichte p und die spezifische Warmekapazitat cq bei konstanter Deformation. Durch 
die Vorgabe von Anfangsbedingungen 


rliz = ro (x), Ur, = 00 (x), Olı=, = Oo (x) 
und einer dazu vertraglichen, disjunkten Vorgabe von Verschiebungen und Spannungen 


r =F (x,t), xel,t>bo 


o'n=ö(x,t), xel.,t2bo 
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auf dem Randgebiet T = T, UT, sowie von Temperatur und Wärmestrom 


0 =Õ (x,t), elo, tr >to 
q:-n=G(x,t), xel,tzbo 
q:n=h (0 (x,t)-Ö(x,t)), xeTog,tzto 


auf dem Randgebiet T = To UT, UTo,, wird das Anfangsrandwertproblem vervollstän- 
digt. n bezeichnet hierbei den nach außen gerichteten Normalenvektor. 

Bei quasi-statischen Lastfällen werden die ursprünglich hyperbolischen Verschiebungs- 
differentialgleichungen (2.10) durch Vernachlässigung der Trägheitsterme pD;v zu 
elliptischen Verschiebungsdifferentialgleichungen, der parabolische Charakter der 
Wärmeleitungsgleichung bleibt erhalten. Für den isotropen Körper lassen sich der 
Wärmeleitfähigkeitstensor K = KI und die Temperaturleitfähigkeit x = K/pca als 
Funktionen von Wärmeleitfähigkeit K, Dichte p und spezifischer Wärmekapazität cq 
bei konstanter Deformation einführen. Es ergeben sich die vollständig gekoppelten 
Verschiebungs- und Wärmeleitungsdifferentialgleichungen 


uAr + (A + u) V(V-1) + pf =(3A +2u) aVO, (2.15) 
(A+2u)e d 90 Q 


Ferenc p (ee) 


KAO 
mit dem dimensionslosen Kopplungsparameter 


_ BAt 2u) a? ToK 
(A + 2u)K 


Das Auflösen von Gleichung (2.14) nach 9 und das Einsetzen in die homogenen 
Gleichungen (2.15) und (2.16) 


3A 4+2 T 
e 


wAr+(A+u+(A+2We)V(V-r)= SS (2.17) 
l+ed 
[a Di 2) s =0 (2.18) 


lässt erkennen, dass die spezifische Entropie s für quasi-statische Problemstellungen 
einer Diffusionsgleichung genügt. Im Falle der einseitigen Kopplung mit € — 0, also 
der Annahme, dass die zeitliche Änderung des Verschiebungsfeldes keinen Einfluss 
auf die Temperatur hat, vereinfacht sich der Lösungsprozess im Allgemeinen deutlich. 
In diesem Fall kann die Wärmeleitungsgleichung unabhängig gelöst werden und das 
erhaltene Temperaturfeld wirkt als Anregung für die Verschiebungsdifferentialglei- 
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chungen. Eine mögliche Begründung für die Vernachlässigung des Kopplungsterms 
liefert beispielsweise die Umformung der homogenen Wärmeleitungsgleichung (2.16) 


zu 
d 
A+2uje SN r 
xA8 = 11+ (A+ 2ue 3 LA 
(3A + 2u) a Ot 
woraus mit den Annahmen EA ri a% < el + 1 und € « 1 eine Vernachlässigung 


des Terms motiviert werden kann. 

Das Lösen der quasi-statischen Wärmeleitungs- und Verschiebungsdifferentialglei- 
chungen, ob nun in der Form (2.15) und (2.16) oder (2.17) und (2.18), ist trotz der 
Linearität oft nur numerisch möglich. Für die Gleichungen wurden verschiedene 
Lösungsmethoden entwickelt, von denen im Folgenden die Helmholtz-Zerlegung und 
die verallgemeinerten Galerkin- und Neuber-Papkovich-Potentiale vorgestellt werden 
sollen. Die Ausführungen orientieren sich an [30, 183-185, 193, 242]. 


e Helmholtz-Zerlegung 

Das Verschiebungsfeld r = rlam + Tel wird als Summe seines lamellaren (wir- 
belfreien) Anteils rjam und solenoidalen (quellenfreien) Anteils r;oı dargestellt. 
Dabei kann der lamellare Anteil riam = V® als Gradient eines Skalars ® und 
der solenoidale Anteil roi = V x W als Rotation eines Vektorfeldes W mit 
V-W = 0 dargestellt werden. Aufgrund V X rlam = 0 und V - rs) = 0 lässt sich das 
Verschiebungsfeld r mit V x r = V X rs) und V-r = V - flam in die beiden Anteile 
aufteilen, sodass sich die Bestimmungsgleichungen V - r = A® und V x r = -AW 
ergeben. Das Einsetzen von 


r=V0+V xY, V-W=0, 
f=Vh+Vxy, V-y=0 


in die Gleichungen (2.15) und (2.16) ergibt 


BA +2uja _ p 
age A+2u ae ory 

p 

Aw — Hy, 

m 

10 (A+2u) cd Q 

A- = . 
| AL EE e K 
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Durch das zusätzliche Eliminieren der Temperatur 0 können die Bestimmungs- 


gleichungen für ® und dh 
Lied _  GBA+2u)aQ P 10 
[a z) Ce (A+2u) K A+2u at an 
Au =- Ca (2.20) 
u 
und durch das zusätzliche Eliminieren von ® können die Bestimmungsgleichun- 
gen für 0 und uh 
Le ës Q p (3A+2p)aTo d 
[a A K A+2u K ane ey) 
Y= gaz (2.22) 


erhalten werden. Bei gegebenen äußeren Anregungen Q und f lassen sich für 
diese Gleichungen zunächst partikuläre Lösungen finden, die Überlagerung einer 
homogenen Lösung sorgt für die Einhaltung der Randbedingungen. Ein Übergang 
zu den einseitig gekoppelten Gleichungen mit € — 0 ist problemlos möglich. 


Verallgemeinerte Galerkin-Potentiale 
Das Verschiebungsfeld r und das Temperaturfeld 9 werden als Kombination eines 
Vektorpotentials F und eines skalaren Potentials Fo 


A+2 2u)e 

DE A Lied AF A+u d A+u+(A+2u)e d vg. F) 
u ot (A + ul Ot 

3A + 2u 

+ aVFo, 

_ (A+2u) ed 


AF 


2 
ag EA H 


(3A + 2u)a x dt 


dargestellt. Durch Einsetzen dieses Ansatzes in die Gleichungen (2.15) und (2.16) 
ergeben sich die vier unabhängigen Gleichungen 


1+eo p 2 

[a <5) aar+ Te (2.23) 
l+eo H Q 

[a S lun + —_ SCH K =0 (2.24) 


für die vier abhängigen Variablen, deren Unabhängigkeit durch die höheren 
Anforderungen an die Differenzierbarkeit erlangt wird. Die erhaltenen Gleichun- 
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gen eignen sich insbesondere für die Berechnung von partikulären Lösungen. Die 
anschließende Überlagerung einer homogenen Lösung sorgt für die Einhaltung 
der Randbedingungen. Der Übergang zu den einseitig gekoppelten Gleichungen 
mit € — 0 ist zwar möglich, hinterlässt aber ein vergleichsweise umfangreiches 
System an Differentialgleichungen. 


Verallgemeinerte Neuber-Papkovich-Potentiale 
Das Verschiebungsfeld r wird als Kombination eines Vektorpotentials und 
eines skalaren Potentials do 


=V(do+x-~)-xX@, Ag =0 (2.25) 


2(A+2u)(1+e) 
A+u+(A+2u)e 
Ansatzes in die homogenen Gleichungen (2.15) und (2.16) ergibt sich 


dargestellt, wobei x = und x = (xyz)' gilt. Durch Einsetzen des 
8 y 8 


Lied 
[a- - 3) Ago = (2.26) 


und die Bestimmungsgleichung für die Temperatur 


_ À+2u 2A +2) Hr od 
LEID (BA +20) a 


Sind die partikulären Lösungen bereits bekannt (beispielsweise durch das Lösen 
der Gleichungen (2.23) und (2.24)) oder keine äußeren Anregungen vorhan- 
den, lassen sich mit diesem Ansatz oft kompakte Lösungen finden. An dieser 
Stelle sei angemerkt, dass der Ansatz (2.25) von Biot [30] ursprünglich für die 
Gleichungen (2.17) und (2.18) verwendet wurde. Durch den Zusammenhang 
zwischen dem Temperatur- und Verschiebungsfeld zur spezifischen Entropie 
gemäß Gleichung (2.14) lässt sich die hier gezeigte Umformulierung aber einfach 
vornehmen (eine entsprechende Herleitung folgt in Unterkapitel 3.1). Der Über- 
gang zu den einseitig gekoppelten Gleichungen mit € — 0 ist in beiden Fällen 
problemlos möglich. Die Neuber-Papkovich-Potentiale sind insbesondere im 
Bereich der vollständig entkoppelten Gleichungen und rein elastischen Lastfällen 
weit verbreitet. 


Der direkte Vergleich zwischen der Helmholtz-Zerlegung und den verallgemeinerten 
Galerkin- und Neuber-Papkovich-Potentialen offenbart zwar eine große Ähnlichkeit der 
Gleichungen (2.19), (2.24) und (2.26), die Interpretation der einzelnen Größen ®, Fo und 
Qo ist allerdings unterschiedlich. In welcher Darstellung die Gleichungen letztendlich 


gelöst werden, hängt von der spezifischen Anwendung ab. 
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2.2.2 Ausgewählte numerische Methoden 


Für das Lösen von partiellen Differentialgleichungen steht heutzutage eine Vielzahl von 
numerischen Methoden zur Verfügung. Erfahrungswerte in der Literatur sprechen im 
Rahmen von kontaktmechanischen Simulationen gegen die Finite-Elemente-Methode 
und für die Rand-Elemente-Methode [175], weshalb Letztere auch in der vorliegenden 
Arbeit zur Anwendung kommt. Im Folgenden werden dafür als wichtig erachtete, 
ausgewählte numerische Methoden aufgeführt und anhand eines Beispiels diskutiert. 
Ausgangspunkt sei zunächst eine partielle Differentialgleichung der Form 


D{r}=p 


mit dem linearen Differentialoperator D{.}, der Inhomogenität p und der gesuchten 
Lösung r. Gelingt das Lösen dieser partiellen Differentialgleichung unter Einhaltung 
der Rand- und Anfangsbedingungen für eine punktförmige Anregung 


Dig} = d(s - š), (2.27) 


so kann die erhaltene Fundamentallösung g anschließend für eine beliebige Anregung p 
durch die kommutative Faltungsoperation 


r(s) = (g * p)(s) = f ecm (2.28) 
R” 


zur gesuchten Lösung r superpositioniert werden. Etwaige Singularitäten von g(s — š) 
werden im Sinne des Cauchy-Hauptwerts behandelt. Das Faltungsprodukt lässt sich 
alternativ als Multiplikation im Bildbereich darstellen, wodurch mit anschließender 
Rücktransformation 


r(s) = 7 {F {g}(@)F {p}(w)} (2.29) 


dasselbe Ergebnis erhalten wird. Hierbei kennzeichnet F{.} die im Allgemeinen 
mehrdimensionale, kontinuierliche Fouriertransformation aus Gleichung (2.1) und 
w bezeichnet die Kreisfrequenzen. Die Transformation von Gleichung (2.27) in den 
Bildbereich bietet zudem häufig einen zweckmäßigen Zugang für das direkte Berechnen 
der Fundamentallösung F {g}(w). 

Die Vor- und Nachteile der Faltungsoperation im Kontext der Kontaktmechanik werden 
beispielhaft anhand der Oberflächenverschiebung eines Halbraums aufgrund eines 
beliebigen Druckfeldes p(s), s = [x, vi auf der Oberfläche erörtert. Eine prinzipielle 
Darstellung des Problems findet sich in Abbildung 2.5(a). Die Verschiebungen an der 
Oberfläche ergeben sich bei gegebenem Druckfeld p(s) aus den Gleichungen (2.28) oder 
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(a) Halbraum mit kontinuierlicher Oberflächenlast (b) Halbraum mit diskretisierter Oberflächenlast 


Abbildung 2.5: Halbraum mit kontinuierlicher und diskretisierter Oberflächenbelastung 
(2.29) entweder mit der bereits 1881 von Boussinesq [119] hergeleiteten Fundamentallö- 
sung im Realraum 


1-v? 1 
nE x2 + y2 


g(s) = (2.30) 


oder alternativ mit der von Ju & Farris [120] und Stanley & Kato [231] berechneten 
Fundamentallösung im Bildbereich 
Lui 2n 


nE j 
wi t wy 


Beide Fundamentallösungen sind im Ursprung singular und für beliebige p(s) muss 


F {g}(w) = 


(2.31) 


in der Regel auf eine numerische Auswertung der Gleichungen (2.28) und (2.29) 
zurückgegriffen werden. Durch die damit verbundene Diskretisierung der beteiligten 
Größen treten zwei Probleme auf. Zum einen werden durch den Übergang zur diskreten 
Fouriertransformation auch ursprünglich nicht-periodische Signale als periodisch 
interpretiert. Eine zunächst denkbare, gemischte Auswertung [144] 


w[s] = DF | {F {g}(w) DF {p}lo]} 


mit der diskreten Fouriertransformation DF {.} ist bei genauerer Betrachtung auf- 
grund des Aliasing-Phänomens in Zusammenhang mit möglicherweise nicht-glatten 
Druckfeldern mit numerischen Fehlern verbunden. Zwar führt eine feinere Diskreti- 
sierung im Frequenzbereich zu einer Verbesserung der Genauigkeit, aber eben auch 
zu überproportional größeren Rechenzeiten sowie Speicherplatzbedarf und ist daher 
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kontraproduktiv. Zum anderen macht die Singularität der Fundamentallösung bei der 
numerischen Auswertung Probleme. Insbesondere das in der Literatur oft verwendete, 
aber mathematisch unsaubere Ersetzen der Singularität lim F {g}(w) entweder durch 
unmittelbare Nachbarwerte oder den lokalen Mittelwert von F {¢}(w) in einem Bereich 
um den Ursprung [150] soll hier vermieden werden. 

Stattdessen wird das ursprünglich kontinuierliche Druckfeld p(s) aus Abbildung 2.5(a) 
durch ein äquidistantes Gitter mit den Elementlangen Ax x Ay als stiickweise konstante 
Funktion P[s] (siehe Abbildung 2.5(b)) diskretisiert. Zusätzlich wird die Fundamental- 
lösung g(s) durch eine Einflussfunktion G(s) ersetzt. Die Einflussfunktion G(s) kann 
durch entsprechende analytische Auswertung von Gleichung (2.28) mit (2.30) und 


a P für — Ax/2 < x < Ax/2 ^ —Ay/2 < y < Ay/2, 
p\s) = 


erhalten werden. Dies liefert den bereits im Jahre 1929 von Love [154] hergeleiteten 


Zusammenhang 
wi, Iw 
=O = E dl sit + Plxax-ax/2 - P|x=x-ax/2 - leen 
Po í Y=y+Ay/2 Y=y-Ay/2 Y=y+Ay/2 Y=y-Ay/2 
=G(s) (2.32) 
mit 


$ = Xin (Y + VX2 +Y?) + Yin (X + VX?+¥2). 


Die Einflussfunktion G(s) enthält keine Singularitäten mehr und das ursprünglich 
kontinuierliche Faltungsintegral in Gleichung (2.28) kann für das nun stückweise 
konstante Druckfeld P[s] als diskrete, lineare Faltung 


+00 +00 


w[s]=(G*P)[s]= 3. 3. GlmAx,nAy|P[(é - m)Ax, (n -n)Ay] (2.33) 


mM=-oo N=—-00 


formuliert werden. Die beschriebene Vorgehensweise verursacht lediglich einen Fehler 
durch die Diskretisierung des Druckfeldes. 
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Gelingt die Umformulierung der diskreten, linearen Faltung in Gleichung (2.33) zu 
einer diskreten, zyklischen Faltung 


M-1N- 
w[s] =(G * P)moals I=), ye IS = M)mod MAX, (n We N)mod nAy]P[mAx, nAy] 


m=0 n=0 


=DF | {DF {G}[w]DF{P}[a]}, (2.34) 


lässt sich das Faltungsprodukt erneut als Multiplikation im Bildbereich ausführen. 
Hierbei bezeichnet mod die Modulo-Operation. Für dienumerische Umsetzung der dis- 
kreten Fouriertransformation kann die schnelle Fouriertransformation genutzt werden. 
Des Weiteren lässt sich das Faltungsprodukt in Gleichung (2.34) auch als Multiplikation 
einer im Allgemeinen vollbesetzten Matrix G und einem Vektor P 


w=GP (2.35) 


darstellen. Für Vektoren w und P der Dimension M ergibt sich für G die Dimension 
M? x M?. Durch Ausnutzen der diskreten, zyklischen Faltung im Bildbereich lässt sich 
also darüber hinaus das speicherplatzintensive Assemblieren der Matrix G vermeiden. 


Da die diskrete Fouriertransformation auf periodischen Funktionen basiert, wird die 
Auswertung von Gleichung (2.33) durch Gleichung (2.34) im Folgenden für ursprüng- 
lich periodische und ursprünglich nicht-periodische Belastungen P [s] anhand eines 
eindimensionalen Beispiels mit P [s] = [Po,...,Pm-ı] und G [s] = [Go,..-,Gu-—a]" 
sowie einer geraden Anzahl M von aquidistanten Elementen mit der Lange As diskutiert. 
Die Vorgehensweise lässt sich ohne Einschränkung auch auf mehrere Dimensionen 
erweitern. 


e P[s] periodisch, G[s] periodisch 
Eine Auswertung von Gleichung (2.34) liefert 


wo Go GM-ı --- Gi Po 
wı G1 Go ... Go Pı 
; = . ; ; . : , (2.36) 
WM-1 Gm-1 Gm-2 --- Go}\Pm-1 


wobei sich die Einträge von G[s] mit Go = G(0), Gi = G(As),..., Guj2 = 
G(MAs/2) und mit Gu-ı = G(-As), Gm-2 = G(-2As), ..., GM/2+1 = G((-M /2+ 
1)As) identifizieren lassen. 

Das Beispiel in Gleichung (2.36) zeigt, dass lediglich die Belastungen von direkt 
angrenzenden Zellen berücksichtigt werden. Es wird daher mit einer Modifikation 
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der Einflussfunktionen gearbeitet. Dafür wird zunächst ein Einflussradius ro; 
definiert, innerhalb dessen der Einfluss eines belasteten Elements berücksichtigt 
wird. Dies ist in Abbildung 2.6 schematisch dargestellt. Für beliebig große 
Einflussradien roi (vgl. Abbildung 2.6(b)) müssen die überlappenden Anteile der 
Einflussfunktion aufsummiert werden und G (s) in Gleichung (2.34) wird durch 
Geff (s) ersetzt. Eine Beschreibung der Auswirkungen eines solchen Einflussradius 
erfolgt in Abschnitt 3.2.4. 


P[s] nicht-periodisch, G[s] periodisch 
Das Rechengebiet wird zunächst künstlich verdoppelt, eine anschließende Aus- 
wertung von Gleichung (2.34) liefert 


wo Go Gom-1 ware Gy Po 
Wy Gi Go ... Go Py 
. ul . . : i f (2.37) 
W2M-1 Gm-ı Gam-2 .-. Go/\Pam-ı 
Für Du Sores PəM-1 = 0, Go = G(0), Gy = G(As), SEA GM = G(MAs) und 


Gom-1 = G(-As), Gom-2 = G(-2As), e.. GM+1 = G((-M + 1)As) sind wo bis 
wm-1 mit dem Ergebnis aus Gleichung (2.33) identisch, wy bis wam-ı können 
nach erfolgter Berechnung verworfen werden. Der beschriebene Algorithmus 
wird bereits vielfach im Rahmen der Kontaktmechanik angewendet [144, 150, 252]. 
Auch für ursprünglich nicht-periodische Funktionen P[s] wird in dieser Arbeit 
mit einem Einflussradius gearbeitet. Für eine korrekte Abbildung des Lastfalls 
ist hierbei ein Einflussradius ro; entsprechend der Größe des Kontaktgebietes zu 
wählen. 


P[s] nicht-periodisch, G[s] nicht-periodisch 

Das Rechengebiet wird zunächst ebenfalls künstlich verdoppelt, eine Auswertung 
von Gleichung (2.34) führt erneut auf Gleichung (2.37). Es werden gleichermaßen 
Py = + = Pam-ı = 0 gewählt, jedoch diesmal Gm = --- = Gam-ı = 0 und 
Go = G(0), Gı = G(As), ..., Gm-ı = G((M - 1)As). Durch diese Vorgehensweise 
wird ein einseitiger Einfluss in positive s-Richtung erreicht. Auch hier können wm 
bis w2m-1 nach erfolgter Berechnung verworfen werden. Diese Variante wird in der 
späteren Arbeit für thermische Probleme von Bedeutung sein, bei denen lediglich 
vergangene Wärmeeinträge einen Einfluss auf das aktuelle Temperaturfeld haben 
dürfen (Kausalität). 


Durch die Nutzung der schnellen Fouriertransformation für die zyklische Faltung 
in Gleichung (2.34) steht sowohl für periodische als auch nicht-periodische Lastfälle 
ein effektives numerisches Werkzeug für die Auswertung von Gleichung (2.33) zur 
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(a) Einflussradius roi, innerhalb dessen der Einfluss (b) Berechnungsschema für Geg(s) bei sich überlap- 
eines belasteten Elements in Gleichung (2.34) penden Einflussfunktionen und beliebig großen 
berücksichtigt wird Einflussradien roi 


Abbildung 2.6: Darstellung periodisch fortgesetzter Oberfläche der Dimension M x N mit verschiedenen 
Einflussradien roi 


Verfügung. Die einzigen Fehlerquellen stellen in allen drei Anwendungsfällen die 
Diskretisierung und die Begrenzung durch den Einflussradius ro; dar. Voraussetzung 
ist, dass die Einflussfunktion G(s) bekannt ist. Die aufgezeigte Vorgehensweise motiviert 
deshalb unter anderem auch die Herleitung von Einflussfunktionen. In der Literatur 
sind diese jedoch lediglich für den rein mechanischen [72], für den rein thermi- 
schen [114] oder für den stationären [138, 161] Lastfall aufgelistet. Einflussfunktionen für 
die voll-gekoppelten, quasi-statischen thermoelastischen Gleichungen und transiente 
Lastfälle wurden nicht gefunden. Es lässt sich an dieser Stelle entsprechend ein weiteres 
Ziel der Arbeit festhalten: 


e Die Herleitung von Einflussfunktionen für rechteckige, thermische und mecha- 
nische Lasten auf der Oberfläche eines Halbraums für die Berechnung des 
Temperaturfeldes und der thermoelastischen Verschiebungen sowie Spannungen 
für die voll-gekoppelten, quasi-statischen Gleichungen des thermoelastischen 
Festkörpers. 


Die bisherigen Ausführungen beschränken sich auf das Vorliegen bekannter Anregun- 
gen p(s), womit sich w(s) durch die diskrete, zyklische Faltung im Bildbereich und 
die Verwendung der schnellen Fouriertransformation zeit- und speicherplatzeffizient 
sowohl für periodische als auch für nicht-periodische Anregungen p(s) berechnen lässt. 
Ist jedoch w(s) bekannt und p(s) von Interesse, müssen inverse Lösungsmethoden 
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angewendet werden. Es wird sich im späteren Verlauf der Arbeit zeigen, dass die Matrix 
G aus Gleichung (2.35) symmetrisch und positiv definit ist und Nebenbedingungen in 
der Form 0 < p < H vorliegen. In der Literatur lassen sich bereits zwei Verfahren finden, 
die denerwähnten Anforderungen entsprechen und die bereits für kontaktmechanische 
Problemstellungen angewendet werden [5, 155, 206, 207, 253]. 

Eine Möglichkeit stellt das Succesive-Over-Relaxation-Verfahren (SOR-Verfahren) aus 
der Familie der Splitting-Verfahren dar, das in Kombination mit einem Moving-Grid 
Ansatz [217] speicherplatzeffizient implementiert werden kann. Eine Konvergenz ist 
bei bekannter Kontaktfläche für positiv definite und symmetrische Matrizen für den 
Überrelaxationsparameter w € (0,2) nachgewiesen. Insbesondere geht das Verfahren 
für w = 1 in das Gauß-Seidel-Verfahren über. Eine Verwendung der diskreten, zykli- 
schen Faltung im Bildbereich lässt das Verfahren jedoch nicht zu. 

Ebenfalls möglich ist das Conjugate-Gradient-Verfahren (CG-Verfahren), für dessen 
Anwendung eine positiv definite und symmetrische Matrix Voraussetzung ist. Es 
konvergiert bei bekannter Kontaktfläche mit monoton fallendem Fehler auf die exakte 
Lösung. Das Verfahren ist auch deshalb interessant, weil es eine Verwendung der 
diskreten, zyklischen Faltung im Bildbereich zulässt. Die Implementierung ist jedoch 
vergleichsweise aufwendiger, da zwei Iterationsschleifen notwendig sind. 

Beide Verfahren können warm gestartet werden. Im Kontext von kontaktmechanischen 
Problemstellungen und großen Gleichungssystemen konvergiert das SOR-Verfahren im 
Vergleich zum CG-Verfahren deutlich langsamer [5, 204, 207, 253]. Eine Implementie- 
rung beider Verfahren hat diese Erfahrungswerte bestätigt, weshalb in Abbildung 2.7 
lediglich ein möglicher Simulationsablauf des CG-Verfahrens für das lineare Gleichungs- 
system Ax = b mit positiv definiter, symmetrischer Matrix A und der Nebenbedingung 
0 < x < H aufgeführt ist. Für eine Implementierung des SOR-Verfahrens wird auf Tian 
und Bhushan [238] verwiesen. 
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Definiere 
Teic = {x € R"|x > OV b> O} 
[pic = 0 


Berechne das Residuum e = b — Ax 
d=e 


Löse Ax = b für x € Teic bis leall et durch 
das Berechnen von: 


Le = er ex/d, Adk 


2. Xk+1 = Xk + ardk 

3. et = Ek — ar Adk 
ER y ag 

4. k= Seat Ek 


5. dky = €k+1 + Brdk 


Berechne g = Ax — b 
Update Teic = {x € R”|x > 0V g <0} 


Tac kon- 
vergiert? 


Update 
T pic = [pic U {x € Teiclx > H} 
Feic = Tac e 


Setze alle x € [pic zu x = H 


nein [pic kon- 


vergiert? 


Abbildung 2.7: CG-Verfahren fiir das lineare Gleichungssystem Ax = b mit Nebenbedingung 0 < x < H 


36 


2.3 Minimalmodelle fiir reibungserregte Schwingungen 


2.3 Minimalmodelle für reibungserregte 
Schwingungen 


Die grundlegenden Ursachen für reibungserregte Schwingungen werden anhand 
zweier Minimalmodelle, die stellvertretend für eine Vielzahl von technischen Systemen 
mit trockener Reibung stehen, beschrieben. Unter reibungserregten Schwingungen 
werden zumeist unerwünschte Schwingungsphänomene verstanden, die auf trockene 
Reibung zurückgeführt werden können, ohne dass dabei weitere äußere Anregungs- 
mechanismen nötig sind. Sie werden deshalb zu den selbsterregten Schwingungen 
gezählt. Die Schwingungen werden häufig eher als Komfort- oder Qualitätsproblem 
anstatt als Ursache für ein Funktionsversagen wahrgenommen. Als Beispiele können 
neben dem Quietschen bei Bremsen oder Straßenbahnen auch das Rupfen bei Kupp- 
lungen oder das Knarzen von Türen genannt werden. Eine entsprechende Vielzahl 
an Forschungsbemühungen lässt sich in der Literatur finden [45, 73, 99-102, 110, 126, 
135, 136, 248]. Im Rahmen dieser Veröffentlichungen haben sich insbesondere zwei 
grundlegende Mechanismen herauskristallisiert, die im Folgenden anhand zweier 
Minimalmodelle diskutiert werden. Die in den Abbildungen 2.8(a) und (b) gezeigten 
Minimalmodelle sind dabei eine abstrakte Vereinfachung der tatsächlichen technischen 
Systeme und erklären jeweils nur mögliche Wege des Stabilitätsverlustes eines sta- 
tionären Systemzustands. Die beiden Instabilitätsmechanismen basieren auf einem 
negativen Reibwertgradienten und einer nicht-konservativen Kopplung, wobei Ersterer 
in Abschnitt 2.3.1 und Letzterer in Abschnitt 2.3.2 diskutiert wird. 


2.3.1 Negativer Reibwertgradient 


Das für die Erklärung niederfrequenter reibungserregter Phänomene wie beispiels- 
weise dem Kupplungsrupfen geeignete Minimalmodell aus Abbildung 2.8(a) besteht 
aus einem Körper mit der Masse m, der durch eine lineare Feder mit der Feder- 
konstanten cı und einen linearen Dämpfer mit der Dämpfungskonstanten dı mit 
der Umgebung verbunden ist. Der Körper reibt auf einem Band, das sich mit der 
konstanten Geschwindigkeit ou bewegt. Die Reibkraft zwischen Band und Körper 
wird durch ein Coulomb’sches Reibgesetz modelliert, der zugehörige Reibwert u(Vre1) 
sei entsprechend den Erläuterungen in Abschnitt 2.1.3 eine Funktion der relativen 
Gleitgeschwindigkeit drei = vo — X. Es ergibt sich die Differentialgleichung für den 
Zustand des Gleitens 


mk + dıX + c1x = U(Vre)Fe SIgN(Yre1), 
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ex 


(a) Minimalmodell fiir negativen Reibwertgradienten (b) Minimalmodell fiir nicht-konservative Kopplung 


Abbildung 2.8: Minimalmodelle fiir reibungserregte Schwingungen 


die sich um die Ruhelage xo = u(vo)Fe/cı linearisieren lässt 


du 
rel 


mAX + |a + ` Ax + c Ax = 0. (2.38) 
Vrel=V0 


Aus Gleichung (2.38) lässt sich bereits erkennen, dass der effektive Dampfungsbei- 


d e $ t d 
5 g Fe für Reibwertgradienten 5 ` - 
Orel Vrel=V0 Orel Vrel=V0 


kann. Eine Auswertung der Eigenwerte 


wert der = dı + 


< 0 negativ werden 


verdeutlicht das Auftreten positiver Realteile für deg < 0. Die so entstehenden rei- 
bungserregten Oszillationen äußern sich meistens als Stick-Slip Schwingungen. Für 
positive Reibwertgradienten wird die Ruhelage des Systems aufgrund von deg > dı 
stabilisiert, was die Bedeutung einer gezielten konstruktiven Beeinflussung des Verlaufs 
der Stribeck-Kurve betont. 

Die Erweiterung des Modells um einen stochastischen Anteil des Reibwerts im Sinne 
eines weißen Rauschens wird für einen einzelnen Freiheitsgrad von Gaus [83] und für 
mehrere Freiheitsgrade von Feng [79] untersucht. 
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2.3.2 Nicht-konservative Kopplung 


Das, für die Erklärung hochfrequenter reibungserregter Phänomene, wie beispielsweise 
dem Quietschen, geeignete Minimalmodell aus Abbildung 2.8(b), besteht aus einem 
Körper mit der Masse m, der durch drei lineare Federn mit den Federkonstanten cy, Cz, 
Ce und zwei linearen Dampfern mit den Dampfungskonstanten dı, da mit der Umgebung 
verbunden ist. Durch die mit dem Winkel y € [0, 7/2] angebrachte lineare Feder mit 
der Federkonstanten c. werden die beiden Freiheitsgrade des Körpers gekoppelt. Der 
Körper reibt an einer durch eine lineare Feder mit der Federkonstanten c, verbundenen 
Kontaktstelle auf einem Band, das sich mit der konstanten Geschwindigkeit vo bewegt. 
Die lineare Feder mit der Federkonstanten c,, kann als Kontaktsteifigkeit interpretiert 
werden [102]. Die Reibkraft wird durch ein Coulomb’sches Reibgesetz modelliert, der 
zugehörige Reibwert u sei konstant. Mit den Beziehungen c11 = Cx + Ce sin? y>®, 
C12 = C21 = Ce sin y cosy > 0 und cz = Cz + Ce cos? y > 0 ergeben sich die gekoppelten 
Differentialgleichungen für den Zustand des Gleitens 


m O)\ {x dı O\/x c11 C12 x\ _ (Ucnz sign(vo — x) 
f CRI A E a Kb (2.39) 


die sich um die Ruhelage 


HCcn=Cc12 
(2) — | cu(C22+Cn)+c21(uCn—c12) E 
= Ge S 


1 
Z 
0 C11 (C22+€n)+C21 (en C12) 


linearisieren lassen 
m 0 Ax di 0 Ax Cii. CI Ucn Ax 0 
= ; 2.40 
D n E D f S be S & C22 + Cp ) be D SS 
4 / 
Aus dem charakteristischen Polynom P(A) = 3. A;A' mit 
i=0 


_ C11 (C22 + Cn) + C21(UCn — C12) 
e 


Ao(u) 72 
A _ dy (C22 + + dacıı SA 
m 
A _(en + C99 + m + didz >0, 
m 

Az fit do > 0, 

m 
Ag =1>0 
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lasst sich nach Liénard und Chipart die Stabilitatsbedingung 
Ay Aı 
A — |A - — 
0 < Ay) < E (A2— 22) 
ableiten, wobei A2— = > Ogilt. Der Stabilitätsverlust entsteht bei Vorzeichenwechsel des 
A 


dritten Hauptminors der Hurwitz-Matrix bei Ao(ui) = = (42 -= 4). Daraus lässt sich 


ableiten, dass der Realteil eines komplex konjugierten Eigenwertpaars das Vorzeichen 
wechselt und eine oszillatorische Instabilität entsteht. 


OI 


w2 


H > > 
Me Mi u He u 
(a) Realteile der Eigenwerte (b) Imaginärteile der Eigenwerte 


Abbildung 2.9: Verlauf von Real- und Imaginärteilen der Eigenwerte für den Sonderfall gleichmafiger 
Dämpfung dı = da 


Für die Anschauung sind die Realteile der komplex konjugierten Eigenwerte für den 
Sonderfall dı = dz als Funktion des Reibwerts u in Abbildung 2.9(a) dargestellt. Sie 
fallen bis zu einem kritischen Reibwert u. zusammen, spalten sich dort auf und 
kurz darauf wechselt einer der beiden Realteile bei u das Vorzeichen, die Ruhelage 
wird instabil. Bei dem kritischen Reibwert u. fallen die in Abbildung 2.9(b) gezeigten 
Eigenkreisfrequenzen zusammen. Dadurch hat sich für diese Art der Instabilität auch 
die Bezeichnung Modenkopplung verbreitet. Es sei jedoch ergänzend angemerkt, dass 
eine Kopplung im strengen Sinne bereits für di # d: nicht mehr vorliegt, es findet 
lediglich eine Annäherung der Realteile und Eigenkreisfrequenzen statt. Des Weiteren 
lässt sich durch den Koeffizienten Ao(u) ableiten, dass sowohl die Kopplung der beiden 
Freiheitsgrade durch die Feder mit der Federkonstanten c. als auch der einseitige 
Einfluss der Reibkraft für das Auftreten dieses Instabilitätsmechanismus notwendig 
ist [101, 102]. 

Für den Übertrag auf reale Anwendungen werden die beiden Freiheitsgrade des Mini- 
malmodells mit zwei Eigenmoden des tatsächlichen technischen Systems identifiziert. 
Der beschriebene Mechanismus lässt sich entsprechend auch in Experimenten beobach- 
ten [213, 249]. Die Änderungen des Reibwerts werden häufig auf verschleißbedingte 
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Veränderungen des Kontakts zurückgeführt. Durch den direkten Vergleich der Mini- 
malmodelle aus Abbildung 2.8 wird außerdem ersichtlich, dass der negative Gradient 
eines geschwindigkeitsabhängigen Reibwerts grundsätzlich auch beim Minimalmodell 
der nicht-konservativen Kopplung zu einer Instabilität führen kann. 

Eine Erweiterung des Modells durch einen stochastischen Anteil des Reibwerts konnte 
nicht in der Literatur gefunden werden. 
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3 Thermomechanischer Kontakt 


gleitender Körper mit rauen 
Oberflächen 


Ein repräsentativer Ausschnitt der scheinbaren Kontaktfläche A, zweier gleitender 
Körper 8; und 82 mit rauen Oberflächen wird in Abbildung 3.1 dargestellt. Die 
Bewegung der beiden Körper ist dabei wie üblich gewissen Randbedingungen auf 
den Randgebieten 08; und 083 unterworfen, die durch die Interaktion mit der Umge- 
bung verursacht werden. Die genauen geometrischen und zeitlichen Eigenschaften 
der tatsächlichen Kontaktfläche A, sind im Allgemeinen zunächst unbekannt. Im 
Folgenden seien nur diejenigen Konfigurationen von Interesse, bei denen eine relative 
Verschiebung der beiden Körper stattfindet, sodass ein gleitender Kontakt mit der 
relativen Geschwindigkeit vrei (X, t) vorliegt. Die Modellbildung konzentriert sich auf 
die Kontaktzone. Das Ziel ist eine möglichst vollumfänglich Abbildung aller Einflüsse 
der dargestellten Konfiguration auf den Reibwert des trockenen Kontaktes, die für die 
Beantwortung der offenen Forschungsfragen notwendig sind. 

Werden die in Kontakt stehenden Körper in der Nähe ihrer Kontaktzone genauer 
untersucht, können dort im Vergleich zum Grundwerkstoff chemisch und mechanisch 
veränderte Mikrostrukturen vorgefunden werden. Letztere lassen sich in plastisch 
verformte und gegebenenfalls verfestigte Bereiche sowie Bereiche mit veränderter 
Gefügestruktur unterteilen [61]. Die feinere Gefügestruktur wirkt sich in der Regel 
auf die mechanischen und thermischen Materialparameter wie beispielsweise Elasti- 
zitätsmodul und Wärmeleitfähigkeit aus, wodurch oberflächennahe Schichten andere 
Werkstoffeigenschaften besitzen als der Grundwerkstoff selbst. Dies ist in Abbildung 3.1 
durch jeweils eine Randschicht angedeutet. 

Für die Berechnung der plastischen Verformung inklusive Verfestigung stehen bereits 
gut entwickelte Methoden zur Verfügung [91, 113, 179, 180], deren zugrundeliegende 
Herleitung auf den Galerkin-Potentialen (vgl. Abschnitt 2.2.1) beruht [53]. Es soll daher 
auf die Abbildung dieses Effektes verzichtet werden, eine nachträgliche Einarbeitung 
kann jedoch erfolgen. Ebenso soll der Einfluss der veränderten Gefügestruktur in 
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Abbildung 3.1: Thermomechanischer Kontakt zweier gleitender Körper Bı und 82 mit rauen Oberflächen — 
ein repräsentativer Ausschnitt der scheinbaren Kontaktfläche A, ist vergrößert dargestellt 


der Nähe der Oberfläche vernachlässigt werden. Entsprechende Veröffentlichungen 
lassen sich ebenfalls in der Literatur finden [195, 196, 205, 216, 225, 226, 235, 254, 
255, 261, 264]. Beide Körper werden in unmittelbarer Umgebung der Kontaktzone 
als homogene, isotrope Halbräume im Rahmen einer linearen, thermoelastischen 
Kontinuumsmechanik modelliert. Das Modell des Halbraums ist im Kontext der 
Kontaktmechanik weit verbreitet [175] und hat seine Anwendbarkeit insbesondere in 
experimentellen Validierungen unter Beweis gestellt [85, 91]. Dennoch sollen zunächst 
die notwendigen Annahmen für die Gültigkeit des Modellierungsansatzes aufgezeigt 
werden: 


1. Es liegen kleine Verzerrungen und kleine Temperaturdifferenzen vor, sodass eine 
lineare Theorie angewandt werden kann. 


2. Die Kontaktfläche ist im Vergleich zur Körpergeometrie und zum Krümmungsra- 
dius der unverformten Oberfläche klein, weswegen der Einfluss der Körpergeo- 
metrie und Körpertemperatur in unmittelbarer Nähe der Kontaktzone vernach- 
lässigbar ist. Dies bedeutet, dass weder Körperform noch Art der Interaktion des 
Körpers mit der Umgebung die Kontaktzone beeinflussen und sich Spannungen 
und Temperaturgradienten auf die Kontaktregion konzentrieren beziehungsweise 
mit größer werdendem Abstand zur Kontaktzone abklingen. 
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3. Die beiden in Kontakt stehenden Körper sind rau, aber nominell glatt. Das 
heißt, dass die jeweilige Amplitude der Rauigkeit im Vergleich zur zugehörigen 
Wellenlänge klein ist. 


In einem homogenen Halbraum können grundsätzlich drei verschieden Wellentypen 
mit unterschiedlichen Wellenausbreitungsgeschwindigkeiten auftreten. Hierzu zählen 


Kompressionswellen vx = „(A + 2u)/p, Scherwellen vs = yu/p und Rayleighwellen 


UR = CUs, wobei ç die Lösung von 
(2—¢*)* = 16(1- ¢?)(1 - c?0¢/v%) 


darstellt und generell vr < vs < vx gilt [119]. 
Weitere Modellierungsannahmen werden wie folgt getroffen: 


4. Das Schwerefeld wird vernachlassigt. 


5. Warme entsteht lediglich durch die Reibleistung, etwaige Warmequellen aufgrund 
von chemischen Reaktionen werden vernachlassigt. 


6. Die relative Gleitgeschwindigkeit ||v,¢1|| der Körper ist im Vergleich zur Aus- 
breitungsgeschwindigkeit der Rayleighwellen vr klein. Dies rechtfertigt die 
Annahme eines quasi-statischen Gleichgewichts und impliziert, dass sich die 
Oberflächenlasten nur langsam mit der Zeit ändern. 


Unter diesen Voraussetzungen können die in Abschnitt 2.2.1 diskutierten quasi- 
statischen Gleichungen für den thermoelastischen Festkörper zur Beschreibung von 
Verschiebungs-, Spannungs- und Temperaturfeldern herangezogen werden. Das Auf- 
bringen instantaner Lasten widerspricht dabei grundsätzlich der Annahme von quasi- 
statischen Lastfällen. Jedoch kann gezeigt werden, dass die resultierende Wellenfront 
schnell genug fortschreitet, sodass bereits nach sehr kurzer Zeit die quasi-statische 
Lösung im Kontaktgebiet gilt [192]. 

Die Oberflächenhöhen im Kontaktbereich sind entweder aus Messungen bekannt oder 
durch Simulationen entsprechend Abschnitt 2.1.2 generiert, wodurch eine diskrete 
Menge an x, y und z-Werten zur Verfügung steht. Die Länge und Breite der Ober- 
flächenelemente sei dabei mit 2a x 2b äquidistant. Es erscheint daher auch aus der 
Sicht der verfügbaren Oberflächendaten analog zu den numerischen Argumenten aus 
Abschnitt 2.2.2 sinnvoll, die voll-gekoppelten und quasi-statischen Gleichungen für den 
thermoelastischen Festkörper zunächst für konstante mechanische oder thermische 
Lasten innerhalb des Gebietes -a < x < a und -b < y < b zu berechnen. Die 
erhaltenen Einflussfunktionen erlauben die Berechnung beliebig geformter und zeitab- 
hängiger mechanischer oder thermischer Oberflächenbelastungen durch eine geeignete 
Superposition (vgl. Abschnitt 2.2.2). Die Berechnung der Einflussfunktionen erfolgt 
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Abbildung 3.2: Ausschnitt eines Halbraums mit äquidistant diskretisierter Oberfläche sowie elementweise 
konstanten mechanischen und thermischen Oberflächenlasten - Druck pa, Scherspannung 
Tox und Toy, Wärmestrom gn und Temperaturdifferenz 94 


in Unterkapitel 3.1 mittels Potentialtheorie. Die erhaltenen Zusammenhänge werden 
daraufhin in aufsteigender Komplexität zunächst in Unterkapitel 3.2 verifiziert, bevor 
die Auswirkungen des Gough-Joule-Effekts in kontaktmechanischen Anwendungen für 
eine Reihe repräsentativer Lastfälle erörtert wird. Der Unterkapitel 3.3 schließt mit einer 
Beschreibung und Diskussion des zu Abbildung 3.1 gehörenden, thermomechanischen 
Modells ab, für das die verschiedenen hergeleiteten Teilmodelle in geeigneter Form 
kombiniert werden. Dabei wird auch der wechselseitige Einfluss zwischen Tempera- 
turerhöhung und Reibleistung berücksichtigt. Die anschließende Formulierung eines 
inversen Problems lässt sowohl bei vorgegebener Bewegungen rı(t) und ro(t) als auch 
bei vorgegebener Normallast eine Berechnung der resultierenden Kontaktgrößen zu. 
Für die Orientierung des Koordinatensystems und die Vorzeichenkonvention sei auf 
Abbildung 3.2 verwiesen. Die Ausführungen der folgenden Unterkapiteln 3.1-3.3 sind 
außerdem bereits in kompakterer Form veröffentlicht [272, 273]. 


3.1 Potentialtheorie 


Die bereits formulierten Annahmen führen auf die homogenen, voll-gekoppelten, quasi- 
statischen Gleichungen (2.15) und (2.16) aus Abschnitt 2.2.1 für den thermoelastischen 
Festkörper. Dementsprechend bieten sich die verallgemeinerten Neuber-Papkovich- 
Potentiale für die weitere Vorgehensweise an. Da Biot diese ursprünglich für die 
Gleichungen (2.17) und (2.18) eingeführt hat [30], wird im Folgenden gezeigt, wie 
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sie sich für die Gleichungen (2.15) und (2.16) einsetzen lassen. Dafür wird der Ansatz 
für die Verschiebungen r in der Form [30] 


r=V(do+x:p)-xp, Ag =0 (3.1) 


mit einem noch unbekannten Parameter x, dem Vektorpotential o. dem skalaren Poten- 
tial do und den kartesischen Raumkoordinaten x = (x y z)" in die voll-gekoppelten 
und quasi-statischen Gleichungen (2.15) und (2.16) eingesetzt 


V ((A + 2u) Ado + (2(A + 2u) - (A + u) x) V - p — (3A + 2u) a0) =0, (3.2) 
A+2u e d _ 06 
nä eer Lider AN 9) SE (3.3) 


Dabei wird ausgenutzt, dass A (x - @) = 2V-@ aufgrund des harmonischen Potentials o 
gilt. Das Weglassen des Gradienten-Operators in Gleichung (3.2) ist gleichbedeutend mit 
der Addition einer Konstanten oder alternativ dem Hinzufügen einer quadratischen 
Funktion zu do oder einer linearen Funktion zu œ. Entsprechendes Auflösen von 
Gleichung (3.2) nach 0 


2 2 2u)- (A 
De Ee, (A +u) Hr. (3.4) 
(3A + 2u) a (BA +2u)a 
und Einsetzen in Gleichung (3.3) ergibt 
An (ru) d 
xAAdo - (1+ €) ER Ado = |2 ER: 2u)“ +(2-yn)e Era P, (3.5) 


wobei der noch freie Parameter x aus der Forderung einer verschwindenden rechten 
Seite von Gleichung (3.5) zu 
_ 2(A+2u)(l+e) 
A+u+(A+2u)e 


bestimmt werden kann. Beim Übergang zu den einseitig gekoppelten Gleichungen mit 
€ — 0 ergeben sich aufgrund von zl. az 4(1 — v) die klassischen Neuber-Papkovich- 
Potentiale. Es sei an dieser Stelle noch ergänzt, dass diese bereits aus der Forderung 
2(A+2u)-(A+ u) x = 0 in Gleichung (3.2) hätten gewonnen werden können. 

Die voll-gekoppelten und quasi-statischen Gleichungen für den thermoelastischen 


Festkörper lassen sich damit durch die Potentiale do = po + IL als Summe der beiden 
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skalaren Potentiale po und II sowie dem Vektorpotential o darstellen, wobei diese die 
folgenden Laplace- und Diffusionsgleichungen 


Ag =0 (3.6) 

Ago =0 (3.7) 
(EE 

AIL se (3.8) 


erfüllen müssen. Der Übergang zu den einseitig gekoppelten Gleichungen kann durch 
€ — 0 erfolgen. Der Spannungstensor ø lässt sich mit 


o=2u [3x (9 gei + VV (ġo +x: p) - [an + (2- ve) d (3.9) 


angeben. Es verbleibt die Aufgabe, genau die Potentialfunktionen o. po und TI zu finden, 
die außer den Gleichungen (3.6)-(3.8) auch die geforderten Randbedingungen erfüllen. 
Dafür werden die Ausdrücke für das Verschiebungs-, Spannungs- und Temperaturfeld 
aus den Gleichungen (3.1), (3.4) und (3.9) zunächst in den Tabellen 3.1 und 3.2 als 
Funktionen der Potentialfunktionen aufgelistet. Der Kopplungsparameter € beeinflusst 
die überwiegende Mehrheit der aufgelisteten Größen. 

Die erhaltenen Lösungen in Tabelle 3.2 eignen sich für Randbedingungen, bei denen 
aus mechanischer Sicht nur eine der drei Spannungskomponenten auf der Ober- 
fläche oa gel, oder o;zz|,-. verbleibt. Darüber hinaus sind die Lösungen 
derart konstruiert, dass die zugehörigen Ausdrücke für die Spannungen keinen 


Kopplungsparameter € enthalten. Aufgrund der voll- Tabelle 3.1: Lösung II 


P P1=M=93=0 


gekoppelten Gleichungen entsteht instantan bei allen drei 


Lastfällen ein Temperaturfeld. Interessant ist der Zusam- po po = 0, #0 
menhang zwischen Spannung o;z|,-.) und Temperaturfeld E oO 

Ou bei Lösung I. Beide sind bis auf die physikalischen R Sti 
Vorfaktoren identisch, besitzen allerdings ein umgekehr- oy 

tes Vorzeichen. Bereits hier lässt sich erkennen, dass im x u 
voll-gekoppelten Fall jede Normalspannungslast mit einer xx 2 GE -AH 
instantanen Temperaturveränderung zusammenhängt, die Oxy 2u ZEN 


im Einklang mit dem Gough-Joule-Effekt bei Druckbelas- oyy 2u E u amt) 
tungen aufgrund der Kompression zu einer Temperaturer- 
höhung und bei Zugbelastungen aufgrund der Expansion 
zu einem Temperaturabfall führt. Für adiabatische Model- Toi SE 
lierungen verbliebe das Temperaturfeld dann konstant, Die Ozz 2u ES = rt) 


Lösung II aus Tabelle 3.1 ist die einzige Möglichkeit, eine o A ATI 


transiente Warmeleitung abzubilden. 
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Nz+ve)laltz+yv)+n+V) 
7 3(z+Y)+ +y) zze OU 0 
zzehe v(nz+vey(a(nz+y)+n+V) EE or u on Ze ala vn 
be HE £ -ZEXE ze _ ZG: nz 229 
voie £ znz dë dbe 
h znz Dye zehe Azo 
T he zzehexe znz =: 2 znz 
cZ 
(Be Gm ir = s 
ce T i ( eZ ` Zeze z) nz Ee znz 
z znz de 550 zech Ah 
be Zezhex ‚he ze (yz + Cz ez) nz D 
EE EE e SE i 
Bee COs : e Ai fezxe (v7 Sieg AA KE xt -ī)- er) nz 70 
E dere (yz _ ) — ee z) niz (£+ be Pe as E Lag | T z Be e 
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Die Lösungen werden im Folgenden durch geeignete Potentialfunktionen für mecha- 
nische Belastungen in Abschnitt 3.1.1 und anschließend für thermische Belastungen in 
Abschnitt 3.1.2 ergänzt und entsprechend den Ausführungen aus Abschnitt 2.2.2 für 
rechteckige Oberflächenlasten innerhalb des Gebietes -a<x<aund-b<y<sbzu 
Einflussfunktionen superpositioniert. Anschließend erfolgt die Superposition für belie- 
bige räumliche oder zeitliche, mechanische und thermische Oberflächenlasten. Dabei 
wird ohne Einschränkung der Allgemeinheit von homogenen Anfangsbedingungen 
der Form 


Thor =0, P| rat =0, rar =0 


ausgegangen. Zudem sollen Verschiebungen, Temperaturen und Spannungen für größer 
werdenden Abstand vom Ursprung gegen null konvergieren 


lim r=0, lim 0=0, lim o=0. 
Ill \|x|| 20 Ile 


3.1.1 Mechanische Oberflächenbelastung 


Die Lösungen I, H und III in Tabelle 3.2 eignen sich für Spannungsrandbedingungen. 
Die Komponenten o, 0, und ozz des Spannungstensors sind unabhängig von €, daher 
können dieselben Potentiale wie im einseitig gekoppelten Fall mit € — 0 verwendet 
werden. Verschiedene Konstruktionsmöglichkeiten für harmonische Funktionen lassen 
sich in [23] finden. Das Potential 


p= “Ta In Lk +y2+22 d 


wird für Lösung I genutzt und erfüllt die Randbedingungen o;x|,-9 = 0, gu 
und 67;|,-9 = -Fzö(x)ö(y), während das Potential 


F 
p SE (zm ( Jt +e +2 +2] dE 


für die Lösungen II und III verwendet wird und die Randbedingungen ozx|,_. = 
—F,6(x)d(y), gel. = -Fyö(x)ö(y) und ozz|,-9 = 0 erfüllt [23]. Beide Potentiale 
genügen Gleichung (3.6). Die physikalischen Vorfaktoren ergeben sich dabei mit r = 
Vx? + y? jeweils aus den Kräftegleichgewichten 


0 


z=0 = 


2rı 


INTE 
0 0 


3.1 Potentialtheorie 


Ozyrdrdy =— Fy, 


ozzrdrdy =-F; 


We E We E 
Ree 


Einige Einflussfunktionen, die eine konstante mechanische Last Tox, Toy oder po 
innerhalb des Gebietes -a < x < a und -b < y < b auf der Halbraumoberfläche mit 
der Temperatur OI A, dem Wärmestrom qz|,_) und der Vertikalverschiebung w|,_9 
auf der Oberfläche in Zusammenhang bringen, werden durch die Auswertung von 
Gleichung (2.28) mit den Lösungen aus Tabelle 3.2 berechnet. In dimensionsloser 
Schreibweise mit x = 2a&, y = 2bn und b = fa können sie als 


8(1 — ul 2L SD as Nese pis (3.10) 
8(1 — ul a -v)e) = -C,,(&,1,B), (3.11) 
szene re =C:z (E,1,B), (3.12) 
re At Temp), 613) 
ssa seat sean ou 
EE on 
Geo 120 das na HEH 6.19) 


geschrieben werden, wobei w = w/2a, Po = a- - v°)po/ nE, Tox = (1 — v?)tox/nE 
und Toy = (1-v 2)tny/nE sowie Czy, ome Cole Tee E T,- und Q, die jeweils 
dimensionslosen Verschiebungen, Spannungen und Einflussfunktionen darstellen. 

Beim Übergang zu den einseitig gekoppelten Gleichungen mit € — 0 verschwindet in 
den Gleichungen (3.13)-(3.16) die Kopplung zur Temperatur und zum Wärmestrom, 
wohingegen die Kopplung zur Vertikalverschiebung in den Gleichungen (3.10) und 
(3.11) für € > œ verschwindet. Des Weiteren geht Gleichung (3.12) für den einseitig 
gekoppelten Fall mit € — 0 in den bereits von Love [154] hergeleitete Zusammenhang 
in Gleichung (2.32) über. Die dimensionslosen Einflussfunktionen sind im Anhang B 
aufgelistet. Weitere Einflussfunktionen für den rein elastischen Lastfall und die einseitig 
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gekoppelten Gleichungen mit € — 0 und a — 0 können Dydo und Busby [72] entnom- 
men werden. Außerdem wird der noch des öfteren benötigte Wert der dimensionslosen 
Einflussfunktion C;; (€, 7, ß) im Ursprung 


Coo = Cx 


so = In( p2+1+ 8) + Bln(/p2 +1+1) 
n=0 


-In(ı/ß? + 1— B) - Bln(a[ß? + 1-1) (3.17) 


hervorgehoben. 


3.1.2 Thermische Oberflachenbelastung 


Die Lösung Il aus Tabelle 3.1 oder Gleichung (3.8) ist äquivalent zur Wärmeleitungs- 
gleichung 


_1+e90 


A : 
S K ot 


(3.18) 


Dies kann durch Einsetzen der nach AIT aufgelösten Gleichung (3.4) in Gleichung (3.5) 
gezeigt werden (unter Berücksichtigung von o = 0). Offensichtlich verkleinert die 
vollständige Kopplung die effektive Temperaturleitfahigkeit, im Folgenden wird Kẹ = 
x/(1+ €) substituiert. 

Gleichung (3.18) soll sowohl für vorgegebene Temperaturfelder als auch für vorgegebene 
Wärmeströme auf der Oberfläche eines Halbraums gelöst werden. In beiden Fällen 
wird die Temperaturverteilung aus der Fundamentallösung der Warmeleitungsglei- 
chung (3.18) für eine instantane Wärmezufuhr g im Vollraum bei x = ž, y = ğ, z = Z 
und t = [47] 


S nz 
S Ok P Ak e(t-t 
SÉ Tee: (3.19) 
8K (mx, (t -t)) 
durch eine entsprechende Superposition im Raum und in der Zeit durch Glei- 
chung (2.28) gefunden. Für die Herleitung der Fundamentallösung sei auf den Anhang C 
verwiesen. Die Fundamentallösung (3.19) lässt bereits erkennen, dass durch die 
kleinere Temperaturleitfähigkeit x, auch stationäre Temperaturfelder vergleichsweise 
später erreicht werden. Bei bekanntem Temperaturfeld kann das Potential II aus 


Gleichung (3.4) 


_1+v 


ATI 
1-v 


ao (3.20) 
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erhalten werden, Verschiebungen und Spannungen folgen aus Tabelle 3.1. 


Für den weiteren Lösungsweg werden zunächst die Ergebnisse von Claesson und Pro- 
bert [58, 59] rekapituliert, die eine Lösung von Gleichung (3.20) für das Temperaturfeld 


VK ¢ 2 
ze ert | 2 Jet Y Jers 
2nKvVt ` \2VKet EIS: 
herleiten, worin erf(...) die Fehlerfunktion bezeichnet. Das Temperaturfeld entsteht 


durch die instantane Wärmequelle ge (x, Y) = qo sign (x) sign (Y) bei ž = 0 und f = 0. 
Die partielle Ableitung dieses Temperaturfeldes 


Ed 


x2 +y2 +22 


d'Ge _ okee Wer 
OxOY 2K (mK et)?!” 


stimmt bis auf den Vorfaktor 1/4 mit der Fundamentallösung aus Gleichung (3.19) für 
ž = ğ = ž = Ťř = 0 überein. Die Autoren erhalten die Lösung von 


x2 +y2 +22 


ES Ie = 1+ v „Ike e e 
drdnu Ly 28 (nxet)? 


durch die allgemeine Lösung von Gleichung (3.20) 


t 
d 1+v 0 Og 
—— ll, = —ak; t+Iht+Ib, 
axdy © an. | Sed ee 


wobei AIT, = 0 gilt [193]. Die beiden Potentiale II, und Ilp werden in diesem Fall so 
gewählt, dass die Spannungen und die Verschiebungen für t > œ erneut gegen null 
konvergieren, womit für die Potentiale 


Il =0, 
1+vaxego 1 


EE [x2 + y2 +z? 


IR = 


und 
Be ELETE 
d Ta 1+vakego 2VKet 


dru ® 1-v nK fx? + y? + 2? 
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folgt. Mit der Definition der Fehlerfunktion 


x 


erf(x) = = J eda 
0 


und der Skalierung s = V4« -t5/./x? + y? + z? lassen sich die verbleibenden Integratio- 
nen in x und y ausfiihren 


1+vaxe.go ` 
Ig = - IZ; = VT wa f 2 ds. 


Spannungen og und Verschiebungen rg für die Wärmequelle ge (X, Y) lassen sich 
wiederum mit Tabelle 3.1 berechnen. 


Zwei verschiedene Superpositionen dieser Lösung werden im Folgenden benutzt. Die 
erste ist eine Superposition zweier Wärmequellen mit der Stärke +qẹ bei z = AZ und 
z=-AZ 
85 = lim (Del. A + Delz=2+12) 
AZ—0 
wat yea) la 
= er erf 
mK Vt Wx et Kies: 


woraus das Potential 


= are 
e EIS? 


; (3.21) 


IL = lim (Helz=2-As + Del Asil 
Az—0 


3222) 


"Sp 


folgt. Erneut können Spannungen og und Verschiebungen rg aus Tabelle 3.1 berechnet 
werden. 
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Die zweite ist eine Superposition einer Warmequelle mit der Stärke +qẹ bei z = AZ und 
einer Wärmequelle mit der Starke -g& bei z = -AZ 


05 = lim (Oelz=z-4ž = Oelz=z+4z) 
AZ—0 
goo 


22 


_ _TK 2802 e Ea, (3.22) 


ert | us Jet Y | 
2 (nxt)? 2VKet 2VK et 


wobei die Größe 9o = qoAz/K bei der Grenzwertbildung konstant gehalten wird. Daraus 
folgt das Potential 


I, = lim (Holz = Del, As 
Az—0 


qo 


1 
_ 1+var.doz | 


2252 
e &e ds. 
1-v Ynk.t 
0 


xs a | ys 
2VKet 2VK et 
Spannungen o} sowie Verschiebungen r, werden auch hier mit Tabelle 3.1 berechnet. 


Vorgegebener Wärmestrom 


Das Temperaturfeld aufgrund einer instantanen Wärmequelle 5 zum Zeitpunkt t = f 
bei x = ž, y = Y auf der Oberfläche eines Halbraums wird durch die Superposition 
zweier instantaner Wärmequellen mit der Stärke +4 bei 2 = AZ und Z = -AZ aus 
Gleichung (3.19) erhalten [47] 


_ @-9°+y-M?+2? 


E IK: e 4xe(t-f) 
AK (nxe (t - f)) 


6, = lim (6 (3.23) 


Až>0 3/2" 


Die Superposition von Gleichung (3.23) für eine transiente, zeitlich begrenzte und 
konstante Wärmequelle qo (t) innerhalb des Gebietes -a < x <aund-b < y < b ergibt 


dl see 4xe(t-F) —| 
AnKVt — 


ST T Oo@|x=x- a — O8|x-= x+a dt 
y-b Y=y+b Y=y-b 


(3.24) 
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mit 


© =e z <=] 
° 2yKevt—t} \oyvecve—F] 


was einer Superposition von Gleichung (3.21) 


entspricht. Folglich führt die Superposition des Potentials IT; 


th 


qa UI) 
ro = H 7 role + | 


xX=x-a 
qo y=y-b y=y+b y 
0 t=t-t t=t-t t=t-t t=t-F 
th Y 
t)1 
09 = ae | ) — Og|x=x+0 + og |x=x-4 — geil — oh |x=x+a dt 
qo 4 +b y=y-b y=y+b =y-b 
0 t=t-f t=t-F t=t-f t=t-f 


zu den gesuchten Verschiebungen und Spannungen. Die Scherspannungen und Verti- 
kalverschiebungen aufgrund einer konstanten, zeitlich begrenzten und rechteckigen 
Wärmequelle qo auf der Oberfläche verschwinden dort, es verbleibt jedoch ein Normal- 
spannungsfeld. Dieses kann in dimensionsloser Form mit x = 2a&, y = 2bn, b = Ba und 
der Péclet-Zahl D. = 4a?/x.At als 


4n?K Oo lz=0 


Zaa(l+v) Ja = Da (&, 1, B, Bt, th) (3.25) 


geschrieben werden, wobei 09,2: = (1- v*)00,z2 [nE und S- die jeweils dimensionslose 
Spannung und Einflussfunktion bezeichnen. Außerdem kann eine dimensionslose Form 
für das Temperaturfeld bei z = 0 aufgrund einer konstanten, zeitlich begrenzten und 
rechteckigen Wärmequelle qn auf der Oberfläche 


2nK OI 


2a qo F Ta (č, N, B; Pe, t, th), (3.26) 
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mit T als dimensionsloser Einflussfunktion hergeleitet werden. 
Die stationären Werte ergeben sich zu 


D 
4Ar?K 0,22 ee 


= lim S =-C 27 
2aa(l+v) qa lim Seal „= Cz (€P), (8.27) 
2nK dä _Iz=0 = 
=lim Ta = 2 
ur qo par q tet Ga (& 1, P) (3 8) 


mit Cz- aus Gleichung (3.12). Eine alternative Herleitungsmöglichkeit bietet das statio- 
näre Potential 


a+ itv &go 
® 1-v4nK 


x2 +y? +z? 


an. Die dimensionslosen Einflussfunktionen können im Anhang B gefunden werden. 
Vollständige Ausdrücke für die stationäre Temperaturverteilung in einem Halbraum 
liefern Manca und Naso [161]. 


Vorgegebenes Temperaturfeld 


Das Temperaturfeld aufgrund einer Oberfläche, die der Temperaturverteilung $ folgt, 
wird durch eine Superpostion einer instantanen Wärmequelle mit der Stärke +5 bei 
Žž = AZ und einer weiteren instantanen Wärmequelle mit der Stärke —ğ bei Z = -AZ aus 
Gleichung (3.19) [47] 


x—¥)2+(y-y)2 +22 
ER ae (3.29) 
8 (nx. (t-9)” 


Il 


erhalten, wobei die Größe S= gAz /K während der Grenzwertbildung konstant gehalten 
wird. Die Superposition von Gleichung (3.29) für ein konstantes, zeitlich begrenztes 
und rechteckiges Temperaturfeld 9, auf der Oberfläche ergibt 


th 


TK ¢ ER | zet H 
ES 3/2 Oelx= Zee + Oelx=x-a - elx= x-a — O@|x=x+a] dt, 
dë (nxe (t - f)) Y=y+ Y=y-b Y=y+b Y=y-b 
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was einer Superposition von Gleichung (3.22) in der Form 


th 


N 2 8 = e 
03 Si 2, Felten + Oglx=x-a — Del - Del 


0 t=t-t 
entspricht. Folglich führt eine entsprechende Superposition des Potentials IT, 


th 


UN REN a shan lal 
ro = — | fo|x=x+a Tal — fo|x=x-a — fo|x=x+a 
Bo 4) Sch ya ` Bach  ® ysy] " 
0 t=t-1 t=t-i t=t-t t=t-f 
th z 
AUTE Ena sl sees lat 
Og = — | Og|x=x+a Og|x=x-a — Og|x=x-a — Og|x=x+a 
Bo Al ya © y=y-b © y=y+b 9 yaycb 
0 t=t-f t=t-t t=t-f t=t-i 


zu den gesuchten Verschiebungen und Spannungen. Die Normalspannungen aufgrund 
einer konstanten, zeitlich begrenzten und rechteckigen Temperaturerhöhung vn auf 
der Oberfläche sind dort null. Es verbleiben jedoch Scherspannungen und Vertikalver- 
schiebungen. Diese werden in dimensionsloser Form mit x = 2aé, y = 2bn, b = Ba und 
der Péclet-Zahl D. = 4a?/x.At als 


An? To,zxlz=0 


aaan Be Ae (3.30) 
Ar? Tezy|,-0 Ea 
Ci NN (3.31) 
4n(1—v) Wol,g sc 
aaa) ore ee) (3.32) 


geschrieben, wobei Tox = (1- v*)o92x/TE, Tozy = (1- v*)oo,2y/ME, Wo = wo/2a 
und Szy sowie S-, und W die jeweils dimensionslosen Spannungen, Verschiebungen 
und Einflussfunktionen bezeichnen. Außerdem kann eine dimensionslose Form für 
das Temperaturfeld bei z = 0 aufgrund eines konstanten, zeitlich begrenzten und 
rechteckigen Temperaturfeldes $y auf der Oberfläche 


ei To (Esp) (8.33) 
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mit T als dimensionsloser Einflussfunktion hergeleitet werden. 
Die stationären Werte ergeben sich zu 


An? GË z=0 _ lim 3 = T (€ ) (3 34) 
alan ee EECH ) 
4n? Tozy =0 T T 
SEA =T ‚Il, H 35 
len 3 lim Szy ve zy (én, B) (3.35) 
_) We = = 
an SY) a eer | See (336) 


all + v) da to that 


mit T;y, e und C,, aus den Gleichungen (3.13), (3.14) und (3.12). Eine alternative 
Herleitungsmöglichkeit bietet sich über das stationäre Potential 


ñ- = 1+vavg Z 
® 1-van Jy +z 


an. Die dimensionslosen Einflussfunktionen sind im Anhang B aufgelistet. 


3.1.3 Superposition für beliebige Oberflächenlasten 


Die erhaltenen Einflussfunktionen aus den Gleichungen (3.10)-(3.16), (3.25)-(3.28) sowie 
(3.30)-(3.36) erlauben die Berechnung beliebig geformter und zeitabhängiger thermi- 
scher oder mechanischer Oberflächenbelastungen durch eine geeignete Superposition. 
Zu diesem Zweck wird ein iteratives Lösungsverfahren für vier häufig verwendete 
Randbedingungen vorgestellt. Während die Fälle 1 und 2 mechanische Belastungen 
der Oberfläche durch Druckfelder abdecken, beinhalten die Fälle 3 und 4 thermische 
Belastungen der Oberfläche im Sinne von Wärmeströmen oder Temperaturfeldern. In 
den letzten beiden Fällen bleibt die Oberfläche spannungsfrei, in den ersten beiden 
Fällen ist die Oberfläche entweder wärmeisoliert oder wärmedurchlässig. 

In allen Fällen ist zunächst eine Diskretisierung der Zeit nötig. Diese erfolgt äquidistant 
mit t = kAt, k = 0,1,...,K und die Verläufe möglicher Wärmeströme qalt) und 
Temperaturen p(t) eines einzelnen Oberflächenelementes sind beispielhaft in den 
Abbildungen 3.3(a) und (b) dargestellt. Da diese nicht nur einen instantanen sondern 
auch zeitlichen Einfluss haben, ist ebenfalls der Verlauf während des Zeitschrittes 
von Bedeutung. Es wird im Folgenden von einem konstanten Verlauf während eines 
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So 


> > 
br tea tk t ben tk-1 tk t 
(a) Wärmestromverlauf (b) Temperaturverlauf 


Abbildung 3.3: Beliebiger zeitlicher Warmestrom- oder Temperaturverlauf eines Oberflachenelementes 


: : s tk tk : : cos : 
Zeitschrittes ausgegangen, sodass für q," und 9," je nach Diskretisierungsvariante 
entweder eine explizite 


q (E n) =), (3.37) 
IE (En) =90 (En, te), (3.38) 
eine implizite 
(En) Et (3.39) 
LE, n) =8y LE, N, tk+1) H (3.40) 


oder eine trapezförmige 


1 

qi" (En) =5 (qo (E;n tk) + qo (EN, Em), (3.41) 
1 

Sa (En) =5 (8n (E, tk) + 80 (E, nern) (3.42) 


Auswertung in Frage kommt. Eine Diskussion zu den Auswirkungen der einzelnen 
Diskretisierungsvarianten erfolgt wahrend der Verifikation in Abschnitt 3.2.2. 


e Fall 1: vorgegebenes Druckfeld, wärmeisolierte Oberfläche 
Die Randbedingungen sind mit 


00 
Ozx|z-9 = 0, |; = 0, Ozz|,<9 = —p(x,y,t), dl =-K ES H =0 
z= 


gegeben und die Vorgehensweise ist in Tabelle 3.3 zusammengefasst. Das vor- 
gegebene Druckfeld p(x, y,t) wird diskretisiert und verursacht in einer ersten 
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Iteration 0 den Temperaturanstieg Op), den Warmestrom qj (p) und das Verschie- 
bungsfeld r(p). Ein Warmestrom q, wird so überlagert, dass qù + q3 = 0 gilt 
und die Normalspannung 07 (q,) und das Temperaturfeld 0(q,) resultieren. Eine 
weitere Normalspannung o7 wird entsprechend of +07 = 0 superpositioniert und 
das Temperaturfeld Oto, ), der Wärmestrom q} (o, und die Verschiebungen r(o7 ) 
folgen. Diese Vorgehensweise wird für n Iterationen wiederholt. Nach jeder 
abgeschlossenen Iteration sind die thermischen Randbedingungen bis auf den 
Warmestrom q} (o,) erfüllt. 


Fall 2: vorgegebenes Druckfeld, wärmedurchlässige Oberfläche 
Die Randbedingungen sind mit 


Ozx|z-0 =0, azyl -o e D. Ozz|,~0 = —p(x, y,t), KI =0 


gegeben und die Vorgehensweise ist in Tabelle 3.4 zusammengefasst. Das vor- 
gegebene Druckfeld p(x, y,t) wird diskretisiert und verursacht in einer ersten 
Iteration 0 den Temperaturanstieg 0} (p) und das Verschiebungsfeld r(p). Ein 
Temperaturfeld 05 wird so überlagert, dass 0; + 09 = 0 gilt und die Scherspannun- 
gen T] (0, ) und die Verschiebungen rt. ) resultieren. Ein Scherspannungsfeld 77 
wird entsprechend € + tT; = 0 superpositioniert und das Temperaturfeld 07 (t7) 
und das Verschiebungsfeld rte, 1 folgen. Diese Vorgehensweise wird für n Ite- 
rationen wiederholt. Nach jeder abgeschlossenen Iteration sind die thermischen 
Randbedingungen bis auf das Temperaturfeld 0; (T,,) erfüllt. 


Fall 3: vorgegebener Wärmestrom, spannungsfreie Oberfläche 
Die Randbedingungen sind mit 


Ozxlz-0 = 0, gel. = 0, Ozz|,-9 = 0, qz|,29 = -K - =q(x,y,t) 
2 |z=0 
gegeben und die Vorgehensweise ist in Tabelle 3.5 zusammengefasst. Der vorge- 
gebene Warmestrom q(x, y,t) wird diskretisiert und verursacht in einer ers- 
ten Iteration 0 den Temperaturanstieg Oto) und die Normalspannung o: (q). 
Eine Normalspannung oŭ wird so überlagert, dass oj + oe = 0 gilt und das 
Temperaturfeld 0(0,), der Wärmestrom q} (07) und die Verschiebungen r(o,) 
resultieren. Ein Wärmestrom qj wird entsprechend qj + ge = 0 superpositioniert 
und die Normalspannung o. (q,) und das Temperaturfeld Ota, 1 folgen. Erneut 
wird eine Normalspannung o, so überlagert, dass of zo = 0 gilt und das 
Temperaturfeld Oo. 1 der Wärmestrom q Lo, ) sowie das Verschiebungsfeld r(o7 ) 
werden erhalten. Diese Vorgehensweise wird fiir n Iterationen wiederholt. Nach 
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jeder abgeschlossenen Iteration sind die mechanischen Randbedingungen bis auf 
ein Normalspannungsfeld o? (q;,) erfüllt. 


e Fall 4: vorgegebenes Temperaturfeld, spannungsfreie Oberfläche 
Die Randbedingungen sind mit 


ge 0A) -9 = 0, Gal = 0, OL =a) 


gegeben und die Vorgehensweise ist in Tabelle 3.6 zusammengefasst. Das vorge- 
gebene Temperaturfeld $(x, y, t) wird diskretisiert und verursacht in einer ersten 
Iteration 0 die Scherspannungen € (9) und das Verschiebungsfeld r(9). Die 
werden so überlagert, dass t$ + T) = 0 gilt und das Tempera- 


0 
turfeld 0j (T7) und die Verschiebungen r(T,) resultieren. Ein Temperaturfeld 0) 


wird entsprechend 6; +0) = Osuperpositioniert und die Scherspannungen 7] (95 ) 
und die Verschiebungen r(0,) folgen. Erneut werden Scherspannungen €, so 
überlagert, dass tT] + t] = 0 gilt und das Temperaturfeld 0f (tĮ) und das 


Scherspannungen T 


Verschiebungsfeld r(t]) werden erhalten. Diese Vorgehensweise wird für n 
Iterationen wiederholt. Nach jeder abgeschlossenen Iteration sind die thermischen 
Randbedingungen mit Ausnahme eines Temperaturfeldes 6} (T,,) erfüllt. 


Für zeitsparende Simulationen können die Iterationen nach Iteration 1 gestoppt werden, 
ansonsten beendet das Erreichen einer vordefinierten relativen oder absoluten Tole- 
ranz cat den Iterationsprozess. Für die einseitig gekoppelten Gleichungen mit € — 0 
ist das Ergebnis bereits nach Iteration 0 für alle Lastfälle exakt. Zudem sind durch 
eine geeignete Superposition der Ergebnisse aus den Abschnitten 3.1.1 und 3.1.2 
auch weitere Randbedingungen umsetzbar. Das vorgestellte Verfahren erlaubt die 
Simulation der quasi-statischen und voll-gekoppelten thermoelastischen Gleichun- 
gen für den Halbraum mit beliebig geformten und zeitabhängigen Lasten auf der 
Oberfläche. Darüber hinaus lassen sich die Ergebnisse leicht auf Anwendungen im 
Rahmen der dreidimensionalen Poroelastizität von Biot [29] übertragen, bei denen die 
zugrundeliegenden Gleichungen strukturell identisch sind. Weitere Informationen zur 
Poroelastizität lassen sich im Anhang A finden. Für detailliertere Informationen zur 
numerischen Umsetzung wird auf die Abschnitte 2.2.2 und 3.2.2 verwiesen. 
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3.1 Potentialtheorie 
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3.22 Verifikation und Diskussion 


Die hergeleiteten Einflussfunktionen werden nachfolgend zunächst in Abschnitt 3.2.1 
verglichen, bevor eine Verifikation mit geeigneten analytischen Lösungen aus der 
Literatur in Abschnitt 3.2.2 erfolgt. Letztere beschränkt sich auf vordefinierte Oberflä- 
chenlasten für die Grenzfalle e > Ooder € > œ, da keine analytischen Lösungen für die 
voll-gekoppelten, quasi-statischen Gleichungen gefunden wurden. Derselbe Abschnitt 
beinhaltet außerdem detaillierte Informationen zu den jeweiligen Superpositionen. 
Bevor die Auswirkungen des Einflussradius im Kontext von periodischen Lasten in 
Abschnitt 3.2.4 erörtert werden, schließt sich eine Diskussion des Einflusses der voll- 
ständigen Kopplung mit € + 0 für verschiedene kontaktmechanische Anwendungsfälle 
in Abschnitt 3.2.3 an. 


3.2.1 Vergleich der Einflussfunktionen 


Der direkte Vergleich der erhalten Einflussfunktionen C-z aus Gleichung (3.12) mit 
T; aus Gleichung (3.26), W aus Gleichung (3.32) und Sz, aus Gleichung (3.25) in 
den Abbildungen 3.4(a)-(d), Tz aus Gleichung (3.13) mit Szy aus Gleichung (3.30) 
in Abbildung 3.4(e) sowie Tay aus Gleichung (3.14) mit oH aus Gleichung (3.31) in 
Abbildung 3.4(f) bietet sich aufgrund der jeweiligen Konvergenz auf die gleichen 
stationären Lösungen an. 

Die Abbildungen 3.4(a), (c) und (d) vergleichen die Einflussfunktionen T;, W und 
S,- mit Cz- zu verschiedenen Zeitpunkten 0 < t < œ entlang der ¢-Richtung bei 
n = 0 aufgrund eines kontinuierlichen Wärmeeintrags qo oder eines kontinuierlichen 
Temperaturfeldes $,. Während Ce T4 und W rein positive Werte besitzen, nimmt oes 
auch negative Werte an. Nach einiger Zeit erreichen alle Kurven eine nahezu stationäre 
Form, deren Werte monoton ansteigen, bis sie für t > © gegen C+- konvergieren. 
Abbildung 3.4(b) zeigt den transienten Verlauf der Einflussfunktionen bei € = 7 = 0. 
Sowohl der Wärmestrom gn als auch die Temperaturerhöhung 9. werden ab dem Zeit- 
punkt fr wieder zu null gesetzt. Offensichtlich folgen Spannung S,- und Verschiebung 
W den Temperaturen T, beziehungsweise T 9, was die Folge ihrer Abhängigkeit vom 
Temperaturanstieg im Halbraum ist. Für t > tp konvergieren Ta, W und Sz wieder 
gegen null. 

Die Abbildung 3.4(e) vergleicht die Einflussfunktion Goz mit den stationären Verläufen 
von (oe und Ts zu verschiedenen Zeitpunkten 0 < t < œ entlang der &-Richtung 
bein = 0 aufgrund eines kontinuierlichen Temperaturfeldes $5. Analog zu den 
Abbildungen 3.4(a)-(d) erreicht die Kurve nach einiger Zeit eine nahezu stationäre 
Form, deren Werte monoton ansteigen, bis sie für t — co gegen T;, konvergiert. 
Die dargestellten Einflussfunktionen enthalten jedoch eine Singularität bei & = 1/2. 
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0 > 
to th tk 
t 
(a) Vergleich der Einflussfunktionen Cz und (b) Zeitlicher Verlauf der Einflussfunktionen (aren 
Ta entlang der €-Richtung bein = Ozu Szz,Tq und W bei € = 7 =0 


verschiedenen Zeitpunkten 


(c) Vergleich der Einflussfunktionen TC; und W ent- (d) Vergleich der Einflussfunktionen Cz, und 
lang der €-Richtung bei ņ = 0 zu verschiedenen Szz entlang der €-Richtung bein = Ozu 
Zeitpunkten verschiedenen Zeitpunkten 


(e) Vergleich der Einflussfunktionen Tz, Czx (f) Vergleich der Einflussfunktionen Leys Gz 
und Szx entlang der €-Richtung bein = 0 zu und $z y entlang der n-Richtung bei € = 0 zu 
verschiedenen Zeitpunkten verschiedenen Zeitpunkten 


Abbildung 3.4: Vergleich verschiedener Einflussfunktionen für B = 1 
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Äquivalente Beobachtungen treffen ebenfalls auf die Verläufe der Einflussfunktionen 
Goa Tay und Cy in Abbildung 3.4(d) entlang der n-Richtung bei € = 0 zu. Da jedoch 
keine Auswertung der Einflussfunktionen auf den Elementgrenzen bei € = n = 1/2 
erfolgt, bereiten die Singularitaten keine Probleme. 


3.2.2 Verifikation 


Die erhaltenen Einflussfunktionen werden mit Ergebnissen von Johnson [119], Jae- 
ger [114], Martini et al. [163], Sternberg und McDowell [232] sowie McNamee und Gib- 
son [167] verifiziert. Die Komplexität der Teilmodelle wird dafür sukzessive erhöht. Für 
die Auswertung der in diesem Kapitel auftretenden diskreten, linearen Faltungen wird 
die diskrete, zyklische Faltung entsprechend den Ausführungen von Abschnitt 2.2.2 
verwendet. Die analytischen Vergleichslösungen aus der Literatur beschränken sich auf 
nicht-periodische, meist radialsymmetrische Oberflachenlasten innerhalb des Radius r, 
sowie auf die beiden Grenzfälle € > 0 und € — oo. Wenn nicht anderweitig angegeben, 
wird im Folgenden eine halbe Elementlänge von a = r,/41 und f = 1 gewählt. 


Mechanisches Teilmodell 


Für ein parabolisches Druckfeld P = Dall r?/ r2)" * innerhalb des Gebietes r < rq 
berechnet Johnson [119] das Verschiebungsfeld 


GN SE (2 — 5) fürr < rfa, ahs 
Po EH ((2-%)arcsi (2) + /5-1 rn 
0 Zi z) arcsin (7 z i 


wobei r den radialen Abstand zum Ursprung bezeichnet. 
Für die einseitig gekoppelten Gleichungen mit € — 0 ergibt eine Superposition im 
Raum von Gleichung (3.12) die dimensionslosen Verschiebungen 


+00 


"la EN > Cz- (m,n,B)P(&E-m,n-n), (3.44) 


m=- nN=—oo 


die aufgrund eines beliebig geformten, dimensionslosen Druckfeldes p(£, n) resultieren. 
Der Vergleich von Gleichung (3.44) mit der analytischen Lösung von Johnson [119] in 
Abbildung 3.5(a) zeigt eine gute Übereinstimmung. Der relative Fehler e,cı zwischen 
den Gleichungen (3.43) und (3.44) ist in Abbildung 3.5(b) dargestellt. 
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A = A f 
120 — Zen -—Johnson 0.25 — Erel in % 
0 


r/Ta r/ta 
(a) Vergleich von Gleichung (3.44) mit Johnson [119] (b) Relativer Fehler e,eı zwischen den Gleichun- 
entlang der r-Richtung gen (3.43) und (3.44) entlang der r-Richtung 


Abbildung 3.5: Verschiebungen aufgrund eines parabolischen Druckfeldes - Vergleich von Gleichung (3.44) 
mit Johnson [119] 


Thermisches Teilmodell - transient 


Für eine rechteckige Wärmequelle go mit der Länge 2l, und der Breite 2l, leitet 
Jaeger [114] das Temperaturfeld 


ve 
v2z2 


gt 
2rıK 8 _ Zäre" =: | 
2a qo d QavV2s 


el xh + Oil, - Oj|x=x-1, - Oylx=r-ı, | ds 
Y=y-I, Y=yt+l, Y=y-lp 


(3.45) 


mit 


e er ll vY 
J = — 
8x25 V8xK2s 


her. Die Warmequelle bewegt sich dabei mit konstanter Geschwindigkeit v in positive 
x-Richtung und Gleichung (3.45) ist im wärmequellefesten Koordinatensystem angege- 
ben. 

Für den beliebigen Wärmestromverlauf eines zunächst einzelnen Elementes in Abbil- 
dung 3.3(a) wird das Temperaturfeld 


2rıK | _ 
2a EA 


tj 
e an Iy ZS 2 Taso (k-i)At (ae SA 
=kA A. ed 


als Superposition in der Zeit von Gleichung (3.26) für die einseitig gekoppelten 
Gleichungen mit € — 0 berechnet. Dafür wird ausgenutzt, dass sich die benötigten 
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ah A — explizit — implizit —trapez 


0 
10 
x / 21, x / 21, 
(a) Vergleich von Gleichung (3.46) mit Jaeger [114] (b) Absoluter Fehler €abs zwischen den Gleichun- 
für die Diskretisierungsvariante (3.41) gen (3.45) und (3.46) für die unterschiedlichen 


Diskretisierungsvarianten (3.37), (3.39) und (3.41) 


Abbildung 3.6: Temperaturfelder aufgrund einer rechteckigen Wärmequelle - Vergleich von Gleichung (3.46) 
mit Jaeger [114] für die verschiedenen Péclet-Zahlen Pe = 0.01,0.1,1 nach einem Gleitweg 
von 6l, 


Einflussfunktionen zu den Zeitpunkten t = kAt für einen Wärmestrom mit der 
Zeitdauer ty, = hAt durch 


"Jl. =kAt " Eule =kAt = 
t„=hAt tn=kAt tn=(k-h)At 


ausdrücken lassen. Die Einflussfunktion Ta ist dadurch in Zeitschritt tų bereits für alle 
vorherigen Zeitschritte bekannt. Die Differenzen Ag = = g = qi, sind identisch null, 
wenn sich der Wärmestrom im jeweiligen Zeitschritt nicht ändert. Die Gesamtlösung 


ergibt sich durch eine zusätzlich Superposition im Raum 


nK +00 +00 k- 
Aa z=0 = > T4 be kAt SE >) T4 t=(k-i)At Aq" zn |- (3.46) 
2a =kA =kA 0 i-1|¢-M 
4 t m=-00 n=— t = H = tn=(k-i)At n-n 
= č=m 
n=n DÉI 


Die Auswertung von Gleichung (3.46) erfordert die Kenntnis der Kontaktflache zu 
jedem Zeitschritt, da dort der Wärmeeintrag stattfindet. 

Gleichung (3.45) wird in Abbildung 3.6(a) nach einem Gleitweg von 61, für die verschie- 
denen Péclet-Zahlen Pe = 0.01,0.1,1 mit Gleichung (3.46) verglichen und zeigt eine 
gute Übereinstimmung. Der Wärmestrom wird dafür entsprechend Gleichung (3.41) 
trapezförmig diskretisiert und die halbe Elementlänge zu a = 1„/41 gewählt. Die 
grauen Bereiche kennzeichnen die Kontaktfläche zu den Zeitpunkten t = to und t = tx. 
In jedem inkrementellen Zeitschritt bewegt sich die Kontaktfläche A; mit A& = 
nach rechts bis sie schließlich ihre aktuelle Position erreicht. Entsprechend wandert 
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auch der Warmeeintrag fiir Gleichung (3.46). Die Kurvenverlaufe zeigen, dass eine 
Temperaturerhohung für höhere Péclet-Zahlen nur dort vorzufinden ist, wo ein Warme- 
eintrag bereits stattgefunden hat. Für kleinere Péclet-Zahlen erhöht sich die Temperatur 
hingegen auch in den anderen Oberflächenbereichen. Der Zusammenhang zwischen 
Gleitgeschwindigkeit und Péclet-Zahl lässt sich somit bereits aus der Anschauung 
erkennen. Im späteren Verlauf der Arbeit beeinflusst die Gleitgeschwindigkeit zudem 
die Menge des Wärmeeintrags, sodass die Absolutwerte der Kurvenverläufe in Abbil- 
dung 3.6(a) dementsprechend zu interpretieren sind. Da außerdem die Kontaktfläche 
Tic (ti) am Ende des jeweiligen Zeitschrittes in der Regel nicht im Voraus bekannt ist, wird 
bei thermomechanischen Kontaktsimulationen ein iteratives Verfahren unvermeidbar 
sein. Die Güte der drei Diskretisierungsvarianten (3.37), (3.39) und (3.41) wird für die 
Peclet-Zahlen P. = 0.01,0.1,1 durch den absoluten Fehler €abs = nK (OI — Ou /ago 
zwischen den Gleichungen (3.45) und (3.46) in Abbildung 3.6(b) verglichen. Der 
Vergleich legt für diesen Bereich an Péclet-Zahlen eine trapezförmige Diskretisierung 
entsprechend Gleichung (3.41) nahe, da hierbei lediglich beim Übergang in und aus 
der Kontaktzone ein Fehler entsteht. 


Thermisches Teilmodell - stationär 


Für eine runde Wärmequelle q (x,y) = go innerhalb des Gebietes r < re leitet John- 


son [119] das stationäre Temperaturfeld 


2nK Ol.-o a ent (3.47) 
SE EEN | 


her, wobei I; und Ip die jeweils kompletten elliptischen Integrale erster und zweiter Art 


bezeichnen. 
Für den Fall einer stehenden Wärmequelle g(&,n) wird das stationäre Temperaturfeld 


+00 


SE a = >. Cy (m,n,B)g(&-m,n-n) (3.48) 


m=- N=—-00 


für die einseitig gekoppelten Gleichungen mit € — 0 als Superposition im Raum von 
Gleichung (3.28) erhalten. 

Die analytische Lösung von Johnson [119] wird in Abbildung 3.7(a) mit Gleichung (3.48) 
verglichen und zeigt eine gute Übereinstimmung. Der relative Fehler e,cı zwischen den 
Gleichungen (3.47) und (3.48) ist in Abbildung 3.7(b) dargestellt. 
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0 — Erel in % 


1 n +> 
0 1 2 3 4 
r/Ta r/Ta 
(a) Vergleich von den Gleichungen (3.48) und (3.56) (b) Relativer Fehler e,eı zwischen den Gleichun- 
mit Johnson [119] entlang der r-Richtung gen (3.47) und (3.48) entlang der r-Richtung 


Abbildung 3.7: Stationäre Temperaturerhöhung aufgrund einer runden Wärmequelle — Vergleich von 
Gleichungen (3.48) mit Johnson [119] 


Thermoelastisches Teilmodell - transient 


Für eine parabolischen Wärmequelle q = qo (1- r?/ r2)" * innerhalb des Gebietes r < r, 
berechnen Martini et al. [163] die thermoelastischen Verschiebungen mittels Faltung 
der zugehörigen Fundamentallösung im Bildbereich numerisch. Die Autoren müssen 
dafür die Diskretisierung im Frequenzbereich verachtfachen, um den numerischen 
Fehler zu minimieren. Der Lastfall entspricht Fall 3 aus Abschnitt 3.1.3 für die einseitig 
gekoppelten Gleichungen mit € — 0. 

Für einen beliebigen Wärmestrom g(&,n,t) werden zunächst die dimensionslosen 
Spannungen 


a 2aa(1+v) S z 
00,2z2| z=0 = > > zz| t=kAt q 
, = 2 tole 
t=kAt An fe ee oe 
n=n 
k-1 
+) Szz| t=(k-i)At Ad: rom | (3-49) 
i=l tn=(k-i)At CS n 
é=m 


n=n 


für die einseitig gekoppelten Gleichungen mit € — 0 als Superposition im Raum und 
in der Zeit von Gleichung (3.25) berechnet. Auch hier wird ausgenutzt, dass sich die 
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_ Zeck  Wolz-o 100 4 Se Jeck  00,z|.-0 
2aall+v) go 2aa(l+v) go 


-— Martini et al. 


15000 
10000 
5000 
0 
0 1 2 3 4 
Wé Wé 
(a) Vergleich von Gleichung (3.50) mit Martini et (b) Zugehörige Spannungen von Gleichung (3.49) 
al. [163] entlang der r-Richtung entlang der r-Richtung 


Abbildung 3.8: Thermoelastische Verschiebungen aufgrund einer parabolischen Wärmequelle - Vergleich 
von Gleichung (3.50) mit Martini et al. [163] für die verschiedenen Zeitdauern des 
Wärmeeintrags tn = 0.05,0.1,0.5,1,5 s 


benötigten Einflussfunktionen zu den Zeitpunkten t = kAt und für einen Wärmestrom 
mit der Zeitdauer tp = hAt durch 


Sux 


t=kAt T Em = Szz t=(k-h)At 
tn=hAt tn=kAt tn=(k-h)At 


ausdrücken lassen. Die Einflussfunktion S,- ist dadurch in Zeitschritt tų bereits für die 
vorherigen Zeitschritte bekannt. Die Differenzen Agi = qi” = Ai sind identisch null, 
falls sich der Wärmeeintrag im jeweiligen Zeitschritt nicht ändert. Die anschließende 
Superposition im Raum von Gleichung (3.12) ergibt die dimensionslosen thermoelasti- 


schen Verschiebungen 


+00 +00 


wol z=0 = >, 2 C-z (m,n,ß) Poza za ëm), (3.50) 


t=kAt a nzað 


die aufgrund eines beliebigen Warmestromverlaufs q(é,17,t) resultieren und für deren 
Berechnung die Kontaktflache in jedem Zeitschritt bekannt sein muss. Als Diskretisie- 
rungsvariante wird auch hier Gleichung (3.41) gewählt. 

Der Vergleich von Gleichung (3.50) mit den numerischen Lösungen von Martini et 
al. [163] in Abbildung 3.8(a) zeigt für die verschiedenen Zeitdauern des Wärmeeintrags 
th = 0.05,0.1,0.5,1,5s und die Temperaturleitfähigkeit « = 3.45e-6 m?/s eine gute 
Übereinstimmung. Die Ergebnisse der Autoren wurden graphisch reproduziert. Die 
zugehörigen Spannungen von Gleichung (3.49) werden im Abbildung 3.8(b) dargestellt. 
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Thermoelastisches Teilmodell - stationär 


Für eine runde Wärmequelle q = go innerhalb des Gebietes r < ra gibt Johnson [119] 
das stationäre, thermoelastische Verschiebungsfeld 


Ark Wo,rel uch -77 (2In (2) +1- 5) fürr < ra, 
2aa(1+v) q ~ Kë er É Gay 
? z KR In (2) 


an, worin rg = (x? + y2)" * einen beliebigen Referenzpunkt auf der Oberfläche bezeich- 
net. 
Die dimensionslosen, stationären thermoelastischen Verschiebungen 


+00 +00 
WO rel Dr = > 2 om (m, n, B) Sa, ZZ E: m = Se, ZZ Ze m (3.52) 


aufgrund einer stehenden Warmequelle g(&,n) ergeben sich für die einseitig gekoppel- 
ten Gleichungen mit € — 0 durch eine Superposition im Raum von Gleichung (3.12). E, 
und n, bezeichnen einen beliebigen Referenzpunkt auf der Oberfläche. Dieser garantiert, 
dass die Summe in Gleichung (3.52) endlich bleibt. Die benötigten dimensionslosen, 
stationären Spannungen 


+00 


A 2 1 
Son] EA S mmBge-mn-n) (3.53) 


m=- N=—oo 


werden durch eine Superposition im Raum von Gleichung (3.27) erhalten. 

Abbildung 3.9(a) vergleicht die analytischen Ergebnisse von Johnson [119] mit Glei- 
chung (3.52) (in diesem Fall wird ro = 4r, gewählt) und zeigt eine gute Ubereinstim- 
mung. Die zugehörigen Spannungen aus Gleichung (3.53) werden in Abbildung 3.7(a) 
dargestellt und fallen mit der stationären Temperaturverteilung aufgrund einer runden 
Wärmequelle zusammen. Der relative Fehler e,eı zwischen den Gleichungen (3.51) und 
(3.52) ist in Abbildung 3.9(b) gezeigt. 

Einen Nachweis für die Divergenz der diskreten, linearen Faltungssumme in Glei- 
chung (3.52) ohne Referenzpunkt liefert beispielsweise das Konvergenzkriterium von 
Bertrand. Dafür wird zunächst eine Abschätzung der stationären thermoelastischen 
Verschiebung aus Gleichung (3.52) ohne Referenzpunkt 


A 


Anik 8) © S > 
== (0,0, 8) = Czz (m,n, B} > Czz (0,0,6) + 
meena OP pies zz (m,n, B} > Caz (0,0, B) Än 
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Wo rel 
_ _4m?K 7e |z=0 T nO 
2aa(1+v) qo 0 Ere] IN 7o 


-— Johnson 


8000 TE — 
4000 
0 —1 > 
0 1 2 3 4 0 1 2 3 
r/Tq r/Tq 
(a) Vergleich von Gleichung (3.52) mit Johnson [119] (b) Relativer Fehler e,eı zwischen den Gleichun- 
entlang der r-Richtung gen (3.51) und (3.52) entlang der r-Richtung 


Abbildung 3.9: Stationäre thermoelastische Verschiebungen aufgrund einer runden Wärmequelle - Vergleich 
von Gleichung (3.52) mit Johnson [119] 


mit ax = 8kC;; (k,k, By getroffen. Mit dem Konvergenzkriterium von Bertrand 
divergiert die Summe >), ax und damit auch die stationäre thermoelastische Verschie- 


bung ol , falls wie im vorliegenden Fall im bk <1mitbx = kIn(k) I — ) —In(k) 


gilt. Dies Get sich durch das Einsetzen E beteiligten Größen leicht zeigen. 


Für ein parabolisches Temperaturfeld 0 = 9o (1 — r2/r2)'/ innerhalb des Gebietes r < ra 
leiten Sternberg und McDowell [232] die stationären thermoelastischen Verschiebungen 


ES ra(1+v) r2 Sé 
wol =E (2-5 fürr <ra, 
= Ge a (3.54) 


v 2 ((2- + Garcia? E 1) 


her. Dieser Lastfall ist vergleichbar mit Fall 4 aus Abschnitt 3.1.3 für die einseitig 
gekoppelten Gleichungen mit € — 0. Mit Ausnahme der physikalischen Vorfaktoren 
fallt das Verschiebungsfeld zudem strukturell mit den Verschiebungen aufgrund eines 
parabolischen Druckfeldes in Gleichung (3.43) zusammen. 
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0.2 — Erel in % 


r/Ta r/Ta 
(a) Vergleich von Gleichung (3.55) mit Sternberg und (b) Relativer Fehler e,eı zwischen den Gleichun- 
McDowell [232] entlang der r-Richtung gen (3.54) und (3.55) entlang der r-Richtung 


Abbildung 3.10: Thermoelastische Verschiebungen aufgrund eines runden Temperaturfeldes - Vergleich von 
Gleichung (3.55) mit Sternberg und McDowell [232] 


Die aufgrund eines stationären Temperaturfeldes 0(&, ul resultierenden, dimensionslo- 
sen, stationären thermoelastischen Verschiebungen 


wo SCH Zare > 2 Czz (m,n,ß)O(&- m,n- n) 
ta 3 > Cax (m, n, B)T GË gah ea) 
1-2 5 
"say 2 3 Czy (m,n, B) Toye] E- mn n) (3.55) 


oo n=- 


werden für die einseitig gekoppelten Gleichungen mit € — 0 durch eine Superposition 
im Raum von den Gleichungen (3.10), (3.11) und (3.36) erhalten. Die benötigten 
dimensionslosen, stationären Scherspannungen 


i _a(1+v) +00 +00 = 

Tone = 2 KS Tex (m,n,B)O(&—m,n—n), (3.56) 
is 1 > 

os| = oc) =”) > $ Fey (mn, 6) OE - m,n—Nn) (3.57) 


resultieren aus einer Superposition im Raum von den Gleichungen (3.34) und (3.35). 

Abbildung 3.10(a) vergleicht die analytischen Ausdrücken von Sternberg und McDo- 
well [232] mit Gleichung (3.55) für die Querkontraktionszahl v = 0.3. Die Ergebnisse 
zeigen eine gute Übereinstimmung. Im Gegensatz zu den stationären thermoelasti- 
schen Verschiebungen aufgrund eines Wärmestroms bleibt das Verschiebungsfeld hier 
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Ur — Erel in % 
—2 
4 > 
0 0.5 1 1.5 2 0 0.5 1 1.5 2 
Ket /r2 Ket [Tq 
(a) Transiente relative Verschiebungen in der Mitte (b) Relativer Fehler e,eı für den Fall einer wärme- 

eines runden Druckfeldes im Vergleich zu durchlassigen Oberflache im Vergleich zu 
McNamee und Gibson [167] Gleichung (3.58) 


Abbildung 3.11: Relative Verschiebungen aufgrund eines parabolischen Druckfeldes fiir den Fall einer 
wärmeisolierten und einer wärmedurchlässigen Oberfläche — Vergleich mit McNamee und 
Gibson [167] 


beschränkt und es wird kein Referenzpunkt benötigt. Der relative Fehler €reı zwischen 
den Gleichungen (3.54) und (3.55) ist in Abbildung 3.10(b) dargestellt. 
Thermoelastisches Teilmodell - transient 


Für den Grenzfall € — oo leiten McNamee und Gibson [167] die transienten, relativen 
Vertikalverschiebungen 


2u 1 AT eal 
org (W = lala = erfe (= +4) = (1 e z) (3.58) 


in der Mitte eines runden Druckfeldes p = po innerhalb des Gebietes r < rą für 


den Fall einer wärmeisolierten Oberfläche und die Querkontraktionszahl v = 0 her. 
Hierbei bezeichnet T die dimensionslose Zeit. Die Autoren berechnen ebenfalls die 
Oberflächenverschiebungen für den Fall einer wärmeisolierten Oberfläche durch nume- 
rische Integration. Sie erhalten die Ergebnisse dabei im Rahmen ihrer Forschungen zu 
fluiddurchströmten, porösen Festkörpern. Die zugrundeliegenden partiellen Differenti- 
algleichungen sind strukturell mit den Verschiebungs- und Wärmeleitungsgleichungen 
identisch, weshalb der hier vorgenommene Übertrag auf thermoelastische Größen 
stattfinden kann. Weitere Informationen hierzu und insbesondere auch zu der von 
McNamee und Gibson [167] angewendeten Theorie des voll-gesättigten, porösen 
Festkörpers können in Anhang A gefunden werden. 

Die vorliegenden Lastfälle entsprechen den Fällen 1 und 2 aus Abschnitt 3.1.3. Für 
die dort beschriebene Vorgehensweise werden die Diskretisierungsvarianten aus 
den Gleichungen (3.41) und (3.42) gewählt und eine halbe Elementlänge von a = 
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r„/101 verwendet. Der Vergleich der erhaltenen transienten, relativen Verschiebungen 
2u(w — pl, al /por. mit den Ergebnissen von McNamee und Gibson [167] in Abbil- 
dung 3.11(a) zeigt für E=1Pa,v = 0, K = 1J/msK, Tp = 1K und a = 1/K eine 
gute Übereinstimmung. Die numerischen Ergebnisse der Autoren wurden graphisch 
reproduziert. Die geringfügigen Unterschiede werden auf die von den Autoren 
erwähnten, numerischen Schwierigkeiten zurückgeführt. Der relative Fehler ereı im 
Vergleich zu Gleichung (3.58) ist für den Fall einer wärmedurchlässigen Oberfläche in 
Abbildung 3.11(b) dargestellt. 


3.2.3 Diskussion des Kopplungseinflusses 


Für die voll-gekoppelten, quasi-statischen Gleichungen konnten keine analytischen 
Lösungen in der Literatur gefunden werden. Da ebenfalls keine numerischen Berech- 
nungen für repräsentative Lastfälle im Kontext der Kontaktmechanik aufzufinden 
waren, soll im Folgenden eine Berechnung durch die in Abschnitt 3.1.3 vorgeschlagene 
Vorgehensweise erfolgen. Dabei ist insbesondere von Interesse, ob eine Vernachlässi- 
gung des Gough-Joule-Effektes durch € — 0 für den hier behandelten Kontakt rauer 
Oberflächen gerechtfertigt ist. 

Eine erste Abschätzung für Lastfälle mit einer Druckbelastung liefert bereits die aus 
Gleichung (3.12) hergeleitete relative Differenz 


"lu Wiig __(1=2V)eE 


= Ku 3.59 
"ml. 14+2(1-v)e (59) 


zwischen den instantanen Verschiebungsfeldern im Falle der einseitigen und der 
voll-gekoppelten Theorie, die bei einer Vernachlässigung des Gough-Joule-Effektes 
mit € — 0 auftritt. Die instantanen Verschiebungen W|;-ọ für eine wärmeisolierte 
und eine wärmedurchlässige Oberfläche sind dabei gleich. Die in Abbildung 3.12 
dargestellten Verläufe zeigen eine Verringerung der instantanen Verschiebung w|,_o 
mit zunehmendem Kopplungsparameter £, wobei die relative Differenz für € — oo 
gegen (1 — 2v) /2 (1 — v) konvergiert. 

Die Materialparameter werden zu E = 210 GPa, v = 0.3, K = 40 J/msK, p = 7850kg/ m’, 
ca = 460 J/kgK und a = 12e-61/K gewählt, woraus der Kopplungsparameter € ~ 
0.011 « 1 mit einer Referenztemperatur Tọ = 293.15 K folgt. Für die aufgeführten 
Materialparameter ergibt sich der relative Fehler aus Gleichung (3.59) zu etwa 0.43 % 
und fällt damit bereits sehr gering aus. Für wärmeisolierte Oberflächen lassen sich 
zudem durch die Gleichungen (3.13)-(3.15) grundsätzliche Abschätzungen der instan- 
tanen Iemperaturerhöhung in Abhängigkeit von der Oberflächenbelastung erhalten. 
Die transienten Auswirkungen können jedoch nur durch numerische Simulationen 
berechnet werden. Durch die verringerte Temperaturleitfähigkeit x. = «/(1+ e) lässt 
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A 2 a mye en 
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Abbildung 3.12: Relative Differenz zwischen den instantanen Verschiebungsfeldern bei voll-gekoppelten und 
einseitig gekoppelten Gleichungen aus Gleichung (3.59) für den Fall einer wärmeisolierten 
und einer wärmedurchlässigen Oberfläche bei Druckbelastung für die Querkontraktions- 
zahlen v = 0,0.1,...,0.4 


sich immerhin schließen, dass die stationären Zustände im Vergleich zu den einseitig 
gekoppelten Gleichungen später erreicht werden. 

Im Folgenden sollen die Fälle 1-4 aus Abschnitt 3.1.3 für nicht-periodische, parabolische 
Lasten vorgestellt werden. Für die numerischen Simulationen wird dafür die halbe 
Elementlänge zu a = r„/101 gewählt, was zu der Péclet-Zahl Pe ~ 0.0392 führt. Erneut 
werden die Diskretisierungsvarianten aus den Gleichungen (3.41) und (3.42) verwendet. 
Als Abbruchkriterium für den iterativen Lösungsprozess wird eine relative Toleranz 
des Verschiebungsfeldes von ci = 10” festgelegt. 

Abbildung 3.13 zeigt einige Ergebnisse für die Fälle 1 und 2 aus Abschnitt 3.1.3 aufgrund 
eines parabolischen Druckfeldesp = po (1 — r?/ r2)" * mit po = 1 GPa sowie ra = 0.5 mm. 
Der zeitliche Verlauf der dimensionslosen Relativverschiebung 2u (w — zl, al /por. in 
der Mitte des Druckfeldes ist in Abbildung 3.13(a) sowohl für eine wärmeisolierte als 
auch für eine wärmedurchlässige Oberfläche dargestellt. Abbildung 3.13(b) zeigt den 
Temperaturanstieg auf der Oberfläche für den Fall einer wärmeisolierten Oberfläche 
entlang der r-Richtung zu den dimensionslosen Zeitpunkten t = 0,0.1,0.5,2. In 
Abbildung 3.13(d) ist die zugehörige transiente Temperaturerhöhung in der Mitte des 
Druckfeldes skizziert, das entstehende Temperaturfeld konvergiert mit fortschreitender 
Zeit gegen null. Im gesamten Halbraum entsteht instantan ein Temperaturfeld, dessen 
Gradient den Verschiebungen entgegenwirkt. An der Oberfläche hat es gemafs dem 
Gough-Joule-Effekt die gleiche Form wie das Druckfeld selbst. Der verbleibende Wär- 
mestrom auf der Oberfläche verändert das Temperaturfeld des Halbraums mit der Zeit 
und verursacht die transienten und zeitverzögerten Verschiebungen. In der Literatur ist 
dieses Phänomen auch als thermoelastische Dämpfung bekannt [33]. Das überlagerte 
Temperaturfeld im Falle einer wärmedurchlässigen Oberfläche verursacht ebenfalls 
zeitverzögerte Verschiebungen, lässt aber die Oberflächentemperatur unverändert. Dies 
führt zu einem geringfügig unterschiedlichen Verschiebungsverlauf. Die dimensionslo- 
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Ket /r2 "ir, 

(a) Transiente relative Verschiebung in der Mitte (b) Temperaturerhöhung entlang der r-Richtung für 
eines parabolischen Druckfeldes für den Fall den Fall einer wärmeisolierten Oberfläche zu den 
einer wärmeisolierten und einer wärmedurchläs- Zeitpunkten 7 = 0,0.1,0.5,2 
sigen Oberfläche 
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(c) Relative Verschiebungen fiir den Fall einer (d) Transiente Temperaturerhöhung in der Mitte 
warmeisolierten und einer warmedurchlassigen eines parabolischen Druckfeldes fiir den Fall 
Oberflache entlang der r-Richtung zu den einer wärmeisolierten Oberfläche 


Zeitpunkten 7 = 0.1,0.5,2 


Abbildung 3.13: Relative Verschiebungen 2u (w - w|;—9) /pora und Temperaturen für die Fälle 1 und 2 aus 
Abschnitt 3.1.3 aufgrund eines parabolischen Druckfeldes 


sen Relativverschiebungen 2u (w — w|,_.) /pora sind für die Zeitpunkte t = 0.1,0.5,2 
und den Fall einer wärmeisolierten und einer wärmedurchlässigen Oberfläche entlang 
der r-Richtung in Abbildung 3.13(c) dargestellt. 

Einige Ergebnisse für die Fälle 3 und 4 aus Abschnitt 3.1.3 aufgrund eines para- 
bolischen Wärmestroms q = qo (1-r?/r2)'” sowie aufgrund eines parabolischen 
Temperaturfeldes 9 = 9o (1- r2/r2)' mit go = 1MJ/sm? und 9o = 100 K sind in 
Abbildung 3.14 dargestellt. Der Temperaturanstieg 9 im Zentrum der Wärmequelle 
wird mit der Temperatur 0|,-ọ für den einseitig gekoppelten Fall in Abbildung 3.14(a) 
verglichen. Die Auswirkungen der verringerten Temperaturleitfähigkeit x. sind deutlich 
sichtbar und führen zunächst zu einer steigenden und dann zu einer sinkenden 
Temperaturdifferenz @|,.., — 0, sobald sich das Temperaturfeld an der Oberfläche 
dem stationären Zustand nähert. Abbildung 3.14(b) stellt die entsprechenden dimen- 
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«t/r2 


(a) Relative Temperaturerhöhung in der Mitte eines (b) Relative Verschiebungen in der Mitte eines 
parabolischen Wärmestroms parabolischen Warmestroms/Temperaturfeldes 


Abbildung 3.14: Relative Temperaturerhöhung und Verschiebungen für die Fälle 3 und 4 aus Abschnitt 3.1.3 


sionslosen Relativverschiebungen K (w — sol. al /@arago im Zentrum der Quelle dar. 
Die Verschiebungen wachsen aufgrund ihrer Abhängigkeit vom Temperaturfeld im 
Halbraum an, das sich im gekoppelten Fall langsamer ausbildet. Des Weiteren sind 
auch die dimensionslosen, relativen Verschiebungen (w — ol. al /og ën für den Fall 
4 in der Abbildung 3.14(b) dargestellt. Da die thermoelastischen Verschiebungen für 
diesen Fall einen stationären Zustand besitzen, kann hier zunächst ein Anwachsen und 
anschließend erneut eine Annäherung der Lösungen beobachtet werden. 

Der Einfluss des Gough-Joule-Effektes auf die Verschiebungs- und Temperaturfelder 
ist in allen vier Beispielen vergleichsweise klein. Der Übergang zu den einseitig 
gekoppelten Gleichungen mit € — 0 erscheint daher für den gewählte Parametersatz 
gerechtfertigt. Für Lastfälle, die eine Druckbelastungen auf der Oberfläche einschlie- 
ßen, kann eine Abschätzung des durch die Vernachlässigung verursachten Fehlers 
mithilfe von Gleichung (3.59) erfolgen. So lassen sich auch für die Werkstoffe Alu- 
minium (7|,-9 — W|ı-o /W|.-9 7 0.89%, € ~ 0.031 bei Tp = 293.15 K) und Kupfer 
(®|.29 — Wrap /Wleag zs 0.59 %, € = 0.019 bei To = 293.15 K) keine nennenswerten 
Unterschiede feststellen. Für transiente Auswirkungen muss auf eine numerische 


Auswertung zurückgegriffen werden. 


3.2.4 Auswirkungen des Einflussradius bei periodischen 
Lastfällen 


Die bisherigen Ausführungen haben sich auf nicht-periodische Lastfälle beschränkt, was 
nicht zuletzt aufgrund von realen Anwendungen sinnvoll erscheint. Bei Modellierungen, 
die einen repräsentativen Ausschnitt der Oberfläche beinhalten, können periodische 
Randbedingungen dennoch von Vorteil sein. Aus diesem Grund sollen im Folgenden die 
Auswirkungen des in Abschnitt 2.2.2 eingeführten Einflussradius roi bei periodischen 
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A —ef — 
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Nr 
(a) Effektive Einflussfunktion rai entlang der (b) Skalierte Anderung Ac T/Coo der effektiven 
€-Richtung bein = 0 fiir verschiedene Einflussfunktion aus Gleichung (3.61) bei € = 7 = 
Einflussradien ro; = NrT mit Nr = 0,5,10,15 0 und größer werdendem Einflussradius roi = 


und einer Periodendauer T = 16 NrT 


Abbildung 3.15: Auswirkung von Periodendauer T und Einflussradius roi = NrT auf die effektive 
—eff 
Einflussfunktion C Se 


— eff 
Lastfällen diskutiert werden. Dafür wird die effektive Einflussfunktion C = für die 
Berechnung der dimensionslosen Oberflächenverschiebungen 


= — eff — 
GU? = om (č, N, b, T, Toi) Pa (3.60) 


verwendet, die die Auswirkungen eines sich mit der Periodendauer T = (Tx, T)" 
wiederholenden rechteckigen Druckfeldes p_, abbildet (vgl. Gleichung (3.12)). Ohne 
Einschränkung der Allgemeinheit wird dabei von einseitig gekoppelten Gleichungen 


mit € — 0 ausgegangen. Die Einflussfunktion ce. hängt zum einen von der dimensi- 
onslosen Periodendauer T und zum anderen vom dimensionslosen Einflussradius ro; 
ab, wobei Letzterer im vorliegenden Beispiel auf ein ganzzahliges Vielfaches der 
Periodendauer Ty = T} = T beschränkt sei. 


Die effektive Einflussfunktion T ist für verschiedene Einflussradien ro; = NrT mit 
Nr = 0,5,10,15 und für die Periodendauer T = 16 entlang der -Richtung bei 7 = 0 
in Abbildung 3.15(a) dargestellt. Die Auswertung erfolgt hierbei nur an den Element- 
mittelpunkten. Es ist erkennbar, dass das Verschiebungsfeld w]|,_, zwar insgesamt für 
einen größer werdenden Einflussradius roi monoton anwächst (und schlussendlich 
gegen unendlich divergiert), die eigentliche Form jedoch nahezu unverändert bleibt. 
Dies wurde bereits in Abschnitt 2.2.2 durch die Singularität der Fundamentallösung 
im Bildbereich in Gleichung (2.31) angedeutet. Dass die Verschiebungen auch für eine 
beliebige Periodendauer T bei größer werdendem Einflussradius roi divergieren, kann 
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anhand des streng monotonen Anwachsens der mit der Periodendauer T skalierten 
Anderung 


—eff —eff — eff 
AC., =C., (E,0,B,T,(n + DT) -C (E,,8,T nT),  n=0,1,... Gei 


bei € = ņ = 0 in Abbildung 3.15(b) gefolgert werden. Gleichung (3.60) lässt außerdem 
darauf schließen, dass für eine konstante Verschiebung des Ursprungs bei einer 


Vergrößerung des Einflussradius ro; ein Druck p_ notwendig ist, der mit 1/ SC skaliert. 
Für ro; > œ konvergiert der Druck py entsprechend gegen null. Der Einflussradius ro; 
hat somit einen direkten Einfluss auf das lokale Druckniveau und damit auf das 
Plastizitätsverhalten. Der in dieser Arbeit eingeführte Einflussradius hat entsprechend 
vergleichbare Auswirkungen wie die beispielsweise von Yastrebov et al. [259, 260] 
vorgeschlagene Wahl von ou > 2r/T für Gleichung (2.6) ohne Verwenden eines 
Einflussradius. 

Da die Einflussfunktionen für thermische Oberflächenbelastungen bei Vorgabe eines 
Wärmestroms gegen Cy konvergieren (vgl. Gleichung (3.28)), konnen die hier getrof- 
fenen Schlussfolgerungen auch auf thermische Problemstellungen tibertragen werden. 
Dementsprechend divergiert das stationäre Temperaturfeld für eine periodische War- 
mequelle qn bei größer werdendem Einflussradius ro; gegen unendlich. 


3.3 Modellbildung und Simulation 


Die bereits in Abschnitt 3.2.2 verifizierten Teilmodelle werden in diesem Unterkapitel für 
den Fall der einseitig gekoppelten Gleichungen mit € — 0 vereinigt. Dies gelingt durch 
die Verbindung der Reibleistung mit dem Wärmeeintrag. Dass eine Vernachlässigung 
des Gough-Joule-Effektes für vergleichbare kontaktmechanische Problemstellungen 
gerechtfertigt ist, haben die Ausführungen in Abschnitt 3.2.3 gezeigt. Die Ergebnisse 
werden dadurch nur geringfügig ungenauer, die Rechenzeitersparnis ist jedoch hoch. 


Ausgangspunkt für das thermomechanische Kontaktmodell sei die in Abbildung 3.16(a) 
definierte Initialkonfiguration. Die beiden Körper mit rauen Oberflächen können sich 
mit ihrer höchsten Erhebung an einer gedachten Ebene berühren. Die Oberflächenhöhen 
zı (x, y) und z2 (x, y) werden in körperfesten Koordinatensystemen beschrieben und 
besitzen die Mittelwerte Sai und Ban, Der Abstand zwischen den Mittelebenen sei 
mit ZH = Z1,max + Z2max gegeben. Jegliche Einflüsse von Temperaturfeldern sind in der 
Initialkonfiguration vernachlässigt. 
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Eine starre, ebene Wand und ein äquivalenter 
rauer Körper mit den Oberflächenhöhen 


z(x,y,t) = zı(x,y) 
+ z2(x = Vrel,xt, H Vrel,yt) (3.62) 


werden in Abbildung 3.16(b) eingeführt. Die 
Höhenwerte der neuen Oberfläche z sind auf- (a) Initialkonfiguration 
grund der relativen Bewegung der beiden glei- 
tenden rauen Körper zeitabhängig. Der Abstand 
zwischen den Körpern in Abbildung 3.16(a) 
und der Abstand zwischen dem äquivalenten 
rauen Körper und der starren, ebenen Wand in 
Abbildung 3.16(b) ist gleich. 

Der äquivalente raue Körper wird durch eine 
Einrückung d in Kontakt mit der starren Wand 
gebracht (siehe Abbildung 3.16(c)). Durch die 
Relativbewegung der beiden Körper und die 
im Kontaktbereich wirkenden Scherspannungen 
kommt es aufgrund von Reibungswärme zu 
den Temperaturerhöhungen 6; im unteren und 


62 im oberen Körper. Infolgedessen treten die 
thermoelastischen Verschiebungen wg = wg, + 
we, auf, die sich aus den einzelnen thermo- (c) Kontaktkonfiguration 
elastischen Verschiebungen der beiden Körper 444,31 dung 2:16: Verschiedene Könfiguralio- 
zusammensetzen. Obwohl der Wärmeeintrag nen für die thermomechani- 
in der gemeinsamen Kontaktfläche stattfindet, sche Kontaktsimulation 

können die Temperaturfelder und thermoelastischen Verschiebungen beider Körper 
dennoch unterschiedlich sein. Die relative Gleitgeschwindigkeit ||vreı|| sei zudem 
klein genug, sodass die Schmelztemperatur nicht erreicht wird. Außerdem wird von 
permanentem Gleiten ausgegangen, wodurch ein partielles Haften ausgeschlossen ist. 


Das überlappende Material 
We(x,y,t)=z+We-ZzH+d (3.63) 


muss durch ein zeitabhängiges, aber zunächst noch unbekanntes Druckfeld p (x, y, t) 
nach unten gedrückt werden. Das in der ebenfalls noch unbekannten Kontaktfläche 
wirkende Druckfeld ist für beide Körper gleich, darüber hinaus wirkende adhäsive 
Anziehungskräfte zwischen den Körpern außerhalb der Kontaktfläche werden ebenso 


82 


3.3 Modellbildung und Simulation 


wie pflügende Verluste vernachlässigt. 
Das lineare, zeitabhängige und dimensionslose Gleichungssystem 


ei 0, (x,y) € Teic, 
| Fz 1 3 Ciz x p + Cut *T + Czy Ty) = Wg u Wp, (x, y) = Joe 
= = 


g, (x, y) € Phic, 
O<p<H, (3.64) 


mit 


1-2v] | Ey (1 = 2v2)(1 + V2) 
Cm 


SE E2 (1 — 2v1)(1 + v1) 


wird aus Abbildung 3.16(c) erhalten. Darin repräsentieren CH *p, Cat, und C, y*Ty die 
diskreten, linearen Faltungen zwischen den Einflussfunktionen und den mechanischen 
Oberflächenlasten. Des Weiteren ergibt sich w; = w,/2a aus Gleichung (3.63) und 
La, Tee und Tnic bezeichnen die Bereiche der Oberflächen, die elastisch, plastisch 
oder nicht in Kontakt stehen. Während des Lösungsprozesses wird der Druck p auf 
positive Werte unterhalb der dimensionslosen Härte H = (1- v?) H/nEı des weicheren 
Materials begrenzt. Das noch überlappende Material wW, = Wy1 + Wp2 im Gebiet Tpic 
wird nach jeder erfolgreichen Lösungsiteration entfernt und volumenerhaltend zu den 
Tälern im Gebiet Phic hinzugefügt [210]. Es muss sichergestellt werden, dass sich die 
Kontaktfläche während dieses Vorgangs nicht verändert. Die gesamte Kontaktfläche 
setzt sich folglich aus Tic = Teic U Tpice zusammen, während eine Lücke g im Gebiet Ta, 
verbleibt. Die Härte H = 3R. lässt sich mit der Streckgrenze Re [238] in Verbindung 
bringen, die ihrerseits eine Funktion der absoluten Temperatur T, der Dehnung € und 
der Dehnrate é ist [157]. Die Streckgrenze R. wird deshalb im Folgenden mittels 


Re(T, €, 8) ~ Ro + Ro0 + RoUrel (3.65) 


approximiert, wobei die beiden Konstanten Rg < 0 und R, > 0 etwaige temperatur- 
und geschwindigkeitsabhängige Effekte berücksichtigten. Entsprechende Messungen 
der Streckgrenze bei unterschiedlichen Temperaturen und Dehnraten lassen sich in [44, 
74, 87] finden. 

Die Wärmeströme q auf den Oberflächen setzen sich aus der Reibungswärme gr und der 
konvektiven Wärmeleitung qn aufgrund der Temperaturdifferenz beider Kontaktkörper 
zusammen. Erstere sind mit der Reibleistung über 


qf = TsUrel (3.66) 
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gekoppelt, wobei ts = R.(T,¢,é)/2 die Versagensscherspannung des weicheren 
Materials bezeichnet. Diese steht über eine lokale Tresca-Vergleichsspannung mit der 
Streckgrenze R. in Zusammenhang. Aufgrund der Abhängigkeit des Wärmestroms q 
vom Temperaturfeld 6 an der Oberfläche wird Gleichung (3.24) zu einer homogenen 
Integralgleichung zweiter Art [132] und unterstreicht die Notwendigkeit eines iterativen 
Lösungsverfahrens. Weiterhin wird angenommen, dass die Reibleistung vollständig in 
Wärme umgewandelt wird [122], wobei eine lineare Wärmeaufteilung schließlich zu 
den einzelnen Wärmeströmen 


_ R - W) ge + h (02 - 61) für Oberfläche 1, a 


y ge+h(01ı - 02) für Oberfläche 2 


führt. Darin bezeichnet y e [0,1] den Wärmeaufteilungskoeffizienten und h den 
Wärmeübergangskoeffizienten. 


Zumindest für ein stationäres Temperaturfeld 6 aufgrund einer stehenden Warme- 
quelle gu = (Ro + Rg@ + R,v)v/2 kann der Einfluss von Gleichung (3.65) analytisch 
bestimmt werden. Eine Verbindung der Gleichungen (3.28) mit (3.65) und (3.66) ergibt 
für die stationäre Temperatur in der Mitte des belasteten Elementes 


A 


boo = e aCoo (Ro + Rov) vV 


=0 = , 
=0 27K — aConpRev 
=0 


woraus das Verhältnis 


wu u, (3.68) 
Ooolo=0 1 - “Sure, 


folgt. Entsprechende Verläufe sind für unterschiedliche Ro und R, in Abbildung 3.17 
dargestellt. Für v — oo konvergiert Gleichung (3.68) auf -2rrKR, /aCooR@Ro. Das 
resultierende Temperaturfeld auf der Oberfläche lässt sich aus Gleichung (3.28) 
berechnen. Außerdem können auch die stationären thermoelastischen Verschiebungen 
hergeleitet werden. Darauf wird hier aber verzichtet. 

Für das Verhältnis Re/Ro folgt ebenfalls die rechte Seite von Gleichung (3.68). Die 
qualitativen Verläufe aus Abbildung 3.17 verdeutlichen somit die Auswirkungen einer 
temperatur- und geschwindigkeitsabhängigen Versagensscherspannung Ts auf die 
Reibkraft. 
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1.5 


0.5 


0 0.25 0.5 0.75 1 
vinm/s 


Abbildung 3.17: Qualitative stationäre Temperaturverläufe aus Gleichung (3.68) bei temperatur- und 
geschwindigkeitsabhängigem Wärmeeintrag für verschiedene Rg und Ry 


Zusammenfassend können die folgenden zusätzlichen thermomechanischen Effekte 
im Vergleich zu rein mechanischen Lastfällen auftreten: 


1. Die thermoelastischen Ausdehnungen aufgrund der Wärmezufuhr an einem 
Asperiten können einen anderen Asperiten in Kontakt bringen. 


2. Die thermoelastisch bedingte Druckerhöhung an einem Asperiten in Kontakt 
kann einen anderen Asperiten den Kontakt verlieren lassen. 


3. Die thermoelastisch bedingte Lasterhöhung kann (weitere) plastische Verformun- 
gen verursachen. 


4. Der Temperaturanstieg aufgrund von Reibungswärme beeinflusst die Materialei- 
genschaften, was wiederum Auswirkungen auf die Kontaktkräfte, die Kontaktflä- 
che und folglich auch den Wärmeeintrag selbst haben kann. 


Die vorgeschlagene Modellierung berücksichtigt somit nicht nur die mechanischen, 
sondern auch die thermischen und thermoelastischen Wechselwirkung aller Asperiten 
von beiden Kontaktkörpern. Die Auswirkung der beidseitigen Rauheit wird durch die 
zeitabhängige äquivalente Oberfläche korrekt abgebildet. 

Für das Gleichungssystem (3.64) sind grundsätzlich zwei verschiedene Lösungsszena- 
rien denkbar. Auf die Vorgabe einer konstanten äußeren Last soll hier verzichtet werden, 
stattdessen wird die Einrückung d vorgegeben. Dies kann mit der Annahme gerecht- 
fertigt werden, dass sich der makroskopische Abstand zwischen den beiden rauen 
Körpern während des Gleitweges über eine Periode der mikroskopischen Rauigkeit 
näherungsweise nicht verändert. Da das Gleichungssystem (3.64) die gleiche Struktur 
wie in rein mechanischen Belastungsfällen besitzt [238], können bereits bewährte 
Lösungsverfahren verwendet werden (vgl. Abschnitt 2.2.2). Abbildung 3.18 fasst den 
gesamten Simulationsablauf zusammen. Die grauen Bereiche heben die im Vergleich 
zu rein mechanischen Simulationen zusätzlich benötigten Simulationsschritte hervor. 
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Materialparameter 
Gemessene Oberflächen © Simulierte Oberflächen 


Berechnung der 
äquivalenten Oberfläche 
aus Gleichung (3.62) 


Berechne 


1. Reı (01, re) und Re (92, 0) nach 
Gleichung (3.65) 


. d für beide Oberflächen nach 


Gleichungen (3.41), (3.66) und (3.67) 
. O1 (q1) und 0z (q2) nach Gleichung (3.46) 
. Woe, (q1) und we, (q2) nach Gleichung (3.50) 


Löse Gleichung (3.64) durch 
CG-Solver aus Abbildung 2.7 
und erhalte Tic, P, Wp1, Wp2, Z 


Tic, 61 und 
02 kon- 
vergiert? 


nein 


ja 


t=t+At Aktualisierung der Oberflachen 
durch w,ı und wy2 


Abbildung 3.18: Simulationsablauf bis gewünschte Zeit teng erreicht ist 
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Kontaktsimulationen 


Um die Auswirkungen der zweiseitigen Rauigkeit auf kontaktmechanische Größen 
herauszuarbeiten, wird das Simulationsmodell aus Kapitel 3 für verschiedene Kombi- 
nationen von Oberflächenstrukturen ausgewertet. Dafür werden zunächst raue Ober- 
flächen mit den Längen Ly = 2mm und Ly = 1mm entsprechend der Vorgehens- 
weise aus Abschnitt 2.1.2 erstellt. Die mit 256 x 128 Elementen äquidistant diskre- 
tisierten Oberflächenhöhen sind normalverteilt Z, ~ N (Sm si und besitzen den 
Mittelwert Sm = Om und die Standardabweichung S, = 0.7 um. Die recht geringe 
Elementanzahl wird gewählt, um nicht nur qualitative, sondern auch quantitative 
Aussagen treffen zu können. Das Leistungsdichtespektrum der erstellten Oberflächen 
folgt der Parametrisierung aus Gleichung (2.6) mit dem Hurst-Exponenten H = 0.8. 
Alle weiteren Oberflächenkennwerte variieren und sind in der Tabelle 4.1 in insgesamt 
sechs Konfigurationen zusammengefasst. Bei den anisotropen Oberflächen wird jeweils 
durch L oder || kenntlich gemacht, ob die zugehörige Oberflächenstruktur transversal 


Tabelle 4.1: Verschiedene Konfigurationen für isotrope und anisotrope Oberflächen 


Konfigurationen wr in mm Wx Wy 

Isotrop 1 Oberflache 1 En 1 1 
Oberfläche 2 Hm 1 1 

Isotrop 2 Oberflache 1 En 1 1 
Oberfläche 2 167 1 1 

Isotrop 3 Oberfläche 1 167 1 1 
Oberfläche 2 167 1 1 

Anisotrop 1 Oberfläche 1 L Yyei 87 1 3 
Oberfläche 2 L Gei 87 1 3 

: Oberfläche 1 || vrer En 3 1 
PON ODS Oberfläche 2 L Gei Ett 1 3 
p Oberfläche 1 || drei En 3 1 
EE Oberfläche 2 Urel En 3 1 
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oder lateral zur Gleitrichtung ausgerichtet ist. 

Die in Tabelle 4.2 aufgefiihrten Materialparameter entsprechen handelstiblichem Stahl 
und werden für beide Kontaktkörper verwendet, um die Untersuchungen auf die Aus- 
wirkungen der zweiseitigen Rauigkeit einzugrenzen. In dem linearen, zeitabhängigen 
Gleichungssystem (3.64) entfallt hierdurch die Kopplung zwischen den tangentialen 
Scherspannungen und den vertikalen Verschiebungen. 


Tabelle 4.2: Verwendete Materialparameter (wenn nicht anderweitig angegeben) 


Parameter Symbol Wert Einheit 
Elastizitätsmodul E 210 GPa 
Querkontraktionszahl v 0.3 - 
Wärmeleitfähigkeit K 40 J/msK 
Dichte p 7850 kg/m? 
Spezifische Warmekapazitat Ca 460 J/kgK 
Wärmeausdehnungskoeffizient a 12e-6 1/K 
Streckgrenze Ro 450 MPa 
Ro 0 MPa/K 
Ry 0 MPas/m 
Warmetibergangskoeffizient h 0 J/Ksm? 
Wärmeaufteilungskoeffizient y 0.5 - 
Referenztemperatur To 293.15 K 


Für die Simulation verschiedener Betriebsbedingungen durchlaufen die Körper 11 
relative Gleitgeschwindigkeiten im Bereich 0 m/s < Vre] < 1 m/s und 21 Einrückungen 
im Bereich zu - 2.5, E + 52, < d < zy. Zudem wird der dimensionslose räumliche 
Einflussradius roi = 113 für alle Einflussfunktionen gewählt. Dies entspricht ungefähr 
0.25 % von Coo aus Gleichung (3.17). Die Berücksichtigung von zeitlichen Einflüssen 
wird auf die vergangenen K = M = 256 Zeitschritte begrenzt. 

Die Darstellung der Ergebnisse erfolgt für die isotropen Oberflächen aus Tabelle 4.1 in 
Unterkapitel 4.1 und für die anisotropen Oberflächen aus Tabelle 4.1 in Unterkapitel 4.2. 
In beiden Unterkapiteln werden zunächst die Auswirkungen eines Einlaufprozesses 
diskutiert, bevor eine schrittweise Auswertung der Simulationsergebnisse stattfindet. 
Für den Einlaufprozess absolvieren alle Oberflächen zunächst einen kompletten Simu- 
lationsdurchlauf mit der Gleitgeschwindigkeit Ga = 1 m/s. Die resultierenden einge- 
laufenen Oberflächen werden anschließend ebenfalls für die Simulationen verwendet. 
Eine Zusammenfassung und vergleichende Gegenüberstellung der Ergebnisse findet 
in Unterkapitel 4.3 statt. Mögliche Korrelationen zwischen den Eigenschaften der 
Oberflächen und dem resultierenden Reibwert werden herausgearbeitet. Abschließend 
behandelt dasselbe Unterkapitel die Erweiterung auf temperaturabhängige Materialpa- 
rameter und deren Auswirkungen. 
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4.1 Isotrope Oberflächen 


Die Verteilungsdichtefunktionen pz(z) der Oberflächenhöhen sind in den Abbildun- 
gen 4.1(a), (c) und (e) für den neuen und eingelaufenen Zustand der isotropen Konfigu- 
rationen aus Tabelle 4.1 dargestellt. In allen drei Fällen sind die Oberflächenhöhen nach 
dem Einlaufprozess linksschief Sr < 0 und sowohl der Wert für die Standardabwei- 
chung S, als auch für die Kurtosis Sr, haben sich verringert. Dennoch verbleiben die 
Oberflächen in den tieferliegenden Bereichen nahezu unverändert und hauptsächlich 
höherliegende Oberflächenbereiche verschleißen. Die beschriebenen Veränderungen 
passen gut zu den Ergebnissen von experimentellen Untersuchungen bezüglich des 
Einlaufverhaltens rauer Körper [116, 124, 198], obwohl die aktuelle Modellierung ein 
vergleichsweise einfaches Verschleißverhalten beinhaltet. Insbesondere das Entstehen 
des charakteristischen Maximums der Verteilungsdichtefunktion wird gut abgebildet. 
Des Weiteren sind die Werte von Kurtosis Sr, und Schiefe S,; in den Abbildungen 4.1(b), 
(d) und (f) für die jeweiligen äquivalenten Oberflächen aus Gleichung (3.62) während 
einer Gleitperiode von Ax = 2 mm sowohl für die neuen als auch für die eingelaufenen 
Oberflächen dargestellt. Die Verläufe zeigen, dass die Oberflächenhöhen der äquivalen- 
ten Oberfläche selbst bei normalverteilten Oberflächenhöhen der Kontaktkörper nicht 
als normalverteilt modelliert werden sollten. Die Erwartungswerte von Standardabwei- 
chung E[S, ], Kurtosis E[Sx.] und Schiefe E[S,;] der jeweiligen äquivalenten Oberflächen 
sind für die Gleitperiode von Ax = 2 mm in Tabelle 4.3 aufgelistet. Sie zeigen, dass die 
äquivalente Oberfläche im Mittel immer noch normalverteilte Höhenwerte besitzt. 

Die Veränderungen in den erhaltenen Verteilungsdichten lassen bereits die Bedeutung 
des Einlaufprozesses für alle nachfolgenden Ergebnisse erahnen. Zudem kann gefolgert 
werden, dass der Einlaufprozess signifikant von den gewählten Betriebsbedingungen 
abhängt. Die neuen und die eingelaufenen Oberflächen werden für die folgenden 
Simulationsergebnisse verwendet. Die Simulationen werden bis zum Erreichen eines 
stationären Betriebszustandes durchgeführt. Die erhaltenen und nachfolgenden Ergeb- 
nisse sind daher periodisch, was auf die periodischen Oberflächen zurückführbar ist. 
Um zunächst einen qualitativen Eindruck der Kontaktgrößenverläufe zu vermitteln, 
werden einige Ergebnisse für die eingelaufenen Oberflächen der Konfiguration Isotrop 1 
für die Gleitgeschwindigkeit v,. = 0.5 m/s und die Einrückung d = 2.08 um in den 
Abbildungen 4.2, 4.3 und 4.4 vorgestellt. Die Abbildungen 4.2 und 4.3 zeigen die 
transienten Temperaturfelder 0; und 62 sowie die transienten thermoelastischen Ver- 
schiebungsfelder woı und wo2 auf den beiden Oberflächen nach einem zurückgelegten 
Gleitweg Ax = 0,0.5,1,1.5 mm in jeweils körperfesten Koordinatensystemen. Die 
Konturlinien für die Temperaturfelder sind dabei in Inkrementen von A@ = 10 K und 
die Konturlinien für die thermoelastischen Verschiebungsfelder in Inkrementen von 
Awg = 0.01 um aufgetragen. Deutlich erkennbar sind die lokalen Oberflächenbereiche 
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A — O1neu und eingelaufen 
—- Oy neu und eingelaufen 


pz(z) 


3 -1.5 0 1.5 
zin um 


(a) Verteilungsdichtefunktionen für Konfiguration 


Isotrop 1 


a "E O; neu und eingelaufen 
—- Oy neu und eingelaufen 


pz(z) 


3 -1.5 0 1.5 
zinum 


(c) Verteilungsdichtefunktionen für Konfiguration 


Isotrop 2 


A — O1neu und eingelaufen 
—- Oy neu und eingelaufen 


pz(z) 


3 -1.5 0 1.5 
z in um 


(e) Verteilungsdichtefunktionen für Konfiguration 


Isotrop 3 


Sku A 
3.5 —-eingelaufen —neu 


3.3 
3.1 + 
a We Bh 
2.7 + 
> 


2.5 
0.7 ES 0.3 01 01 0.3 
Ssk 


(b) Kurtosis Sku und Schiefe Ssg für Konfiguration 
Isotrop 1 


Sku A 
3.5 —-eingelaufen —neu 


3.3 7 
3.1 d 
2.9 
2.7 A 


2.5 + + + | >> 
0.7 -0.5 -03 -01 01 03 


Ssk 


(d) Kurtosis Sku und Schiefe Ssg für Konfiguration 
Isotrop 2 


Sku 

3.5 -- en —neu 
3.3 

3.1 

2.9 7 

2.7 


2.5 — > 
0.7 -0.5 -0.3 -0.1 01 03 


Ssk 


(f) Kurtosis Sku und Schiefe Ssp für Konfiguration 
Isotrop 3 


Abbildung 4.1: Ausgewählte Eigenschaften der neuen und eingelaufenen Oberflächenhöhen von den 
Konfigurationen Isotrop 1, Isotrop 2 und Isotrop 3 aus Tabelle 4.1 
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mit starkem Temperaturanstieg, die sich mit dem Gleitweg ändern und durch den 
kurzzeitigen Kontakt tibereinandergleitender Asperiten entstehen. Die Konturlinien 
fiir die thermoelastischen Verschiebungen zeigen kleinere Gradienten als die Kon- 
turlinien für die Temperaturfelder, was auf die Abhängigkeit der thermoelastischen 
Verschiebungen von den Temperaturfeldern innerhalb der Körper zurückzuführen 
ist. Die transienten Temperatur- und thermoelastischen Verschiebungsverläufe an der 
Position x = 1.25 mm und y = 0.75 mm sind in Abbildung 4.4(a) während des gesamten 
Gleitwegs von Ax = 2 mm dargestellt. Dieser Punkt ist in den Abbildungen 4.2 und 4.3 
jeweils mit einem roten Quadrat markiert und bewegt sich auf der Oberfläche 2 aufgrund 
der Relativbewegung der Kontaktkörper entlang der x-Richtung. Die grauen Bereiche 
in Abbildung 4.4(a) kennzeichnen die Zeiträume, zu denen die beiden Oberflächen 
an diesem Punkt in Kontakt stehen. Auch hier lassen sich erneut die auftretenden, 
kurzzeitigen Temperaturerhöhungen, sogenannte Blitztemperaturen (vgl. Abschnitt 1.2), 
deutlich erkennen. Die dargestellten Verläufe hängen einerseits von der Kontakthistorie 
an diesem Punkt, aber auch vom Kontaktverhalten der umgebenden Asperiten ab. 

In Abbildung 4.4(b) sind die Verläufe der Temperatur 04 = %0/MN, des Nenn- 
druckes pa = I p/MN, der Scherspannung TA = 3. Ts/MN und des daraus resultie- 
renden Reibwerts ua = "alba erkennbar. Die Werte oszillieren mit fortschreitendem 
Gleitweg aufgrund der beidseitigen Rauigkeit und der konstanten Einrückung. Die 
Erwartungswerte dieser Kurven werden später für verschiedene Gleitgeschwindig- 
keiten und Einrückungen miteinander verglichen. Die gestrichelten Linien zeigen 
darüber hinaus die gleichen Kontaktgrößen für den Fall einer temperaturabhängigen 
Streckgrenze Ro = -1 MPa/K (vgl. [44, 74, 87]). Der Druck pa sinkt aufgrund der 
geringeren Reibleistung und den daraus resultierenden kleineren Temperaturen und 
thermoelastischen Verschiebungen, während die Scherspannung t, aufgrund der 
sinkenden Versagensscherspannung reduziert wird. Dies führt im vorliegenden Fall zu 
einem kleineren Reibwert. Die gepunkteten Linien repräsentieren überdies die gleichen 
Kontaktgrößen für den Fall einer geschwindigkeitsabhängigen Streckgrenze R, = 
100 MPas/m. Erst ab Werten dieser Größenordnung treten vergleichbare Auswirkungen 
wie bei einer temperaturabhängigen Streckgrenze auf. Der Druck p4 steigt aufgrund 
der höheren Reibleistung und den daraus resultierenden größeren Temperaturen und 
thermoelastischen Verschiebungen, während die Scherspannung t, aufgrund der 
wachsenden Versagensscherspannung erhöht wird. Dies führt im vorliegenden Fall 
zu einem größeren Reibwert. Werden die Effekte von temperatur- und geschwindig- 
keitsabhängiger Streckgrenze mit Rg = -1 MPa/K und Ry = 100 MPas/m überlagert, 
ergibt sich größtenteils eine Auslöschung der Auswirkungen, weshalb die erhaltenen 
Kurven aus Gründen der Übersichtlichkeit nicht eingezeichnet sind. Zudem gestaltet 
sich eine Abschätzung von R, aus Messdaten recht schwierig, da diese im Allgemeinen 
in Abhängigkeit der Dehnrate gegeben sind (vgl. [44, 74, 87]). Deshalb erfolgt bei den 
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= Se 


GA 


xinmm 
(a) mol. a nach einem Gleitweg Ax = 0 mm 


(b) 6|,-9 nach einem Gleitweg Ax = 0 mm 


—, Awe = 0.01 um >” 


nmm 


um 


ee 


xinmm 


(c) mol. a nach einem Gleitweg Ax = 0.5 mm 


(g) vol. a nach einem Gleitweg Ax = 1.5mm (h) OI a nach einem Gleitweg Ax = 1.5mm 


Abbildung 4.2: Temperaturfelder O und thermoelastische Verschiebungsfelder wg auf Oberfläche 1 für 
unterschiedliche Gleitwege 
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1 — Awọ = 0.01 um 
\O 2 


x in mm 


x in mm x in mm 
(c) OI. a nach einem Gleitweg Ax = 0.5 mm (d) pel. a nach einem Gleitweg Ax = 0.5 mm 


— Awọ = 0.01 um 
I (o, 


x in mm 


(f) mol. a nach einem Gleitweg Ax = 1 mm 
— Awọ = 0.01 um 


14 LOSE 


— 49 =10K 


1.5 


xinmm 


(g) OI. a nach einem Gleitweg Ax = 1.5 mm (h) ol. a nach einem Gleitweg Ax = 1.5mm 


Abbildung 4.3: Temperaturfelder O und thermoelastische Verschiebungsfelder wg auf Oberfläche 2 für 


unterschiedliche Gleitwege 
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meisten Ergebnissen eine Beschränkung auf temperaturabhängige Effekte. 

Die Entwicklung des Temperaturfeldes 0 bei gleich bleibender Einrückung und unter- 
schiedlichen Gleitgeschwindigkeiten lässt sich gut anhand der Verteilungsdichtefunk- 
tion pe(0) beobachten. Diese ist für die Gleitgeschwindigkeiten vre, = 0.1,0.5,1m/s 
in Abbildung 4.4(d) zu sehen. Durch die auftretenden Blitztemperaturen entstehen die 
Ausläufer der Verteilungsdichtefunktionen hin zu höheren Temperaturen. Außerdem 
lässt sich erkennen, dass die Varianz der auftretenden Temperaturwerte mit steigender 
Gleitgeschwindigkeit zunimmt. 

Die Auswirkungen eines konvektiven Wärmeaustausches zwischen den Kontaktkör- 
pern in der tatsächlichen Kontaktfläche auf die Temperatur und die thermoelastischen 
Verschiebungen sind in Abbildung 4.4(c) dargestellt. Hierfür werden die Simulationen 
mit dem Wärmeübergangskoeffizienten h = 1 MJ/Ksm? wiederholt und die absoluten 
Änderungen in den erhaltenen Temperatur- und thermoelastischen Verschiebungsver- 
läufen wiedergegeben. Die Auswirkungen lassen sich erst bei diesem vergleichsweise 
hohen Wärmeübergangskoeffizient beobachten und fallen recht gering aus. Für h — oo 
sind die lokalen Temperaturen der Körper in der tatsächlichen Kontaktfläche zu 
jedem Zeitpunkt identisch. Der Wärmeübergangskoeffizient bietet damit eine gute 
Möglichkeit, die Simulationsergebnisse an Experimente anzupassen. 

In den Abbildungen 4.4(e) und (f) sind die Bereiche der beiden Oberflächen in schwarz 
markiert, die während des gesamten Gleitwegs von Ax = 2 mm mindestens einmal 
in Kontakt stehen. Bei beiden Oberflächen sind dies in etwa 31 %, wohingegen die 
tatsächliche Kontaktfläche zu jedem inkrementellen Zeitpunkt nur etwa 6 % beträgt. 
Somit sind deutlich mehr Oberflächenbereiche am eigentlichen Reibprozess beteiligt, 
als rein statische Simulationen vermuten lassen würden. 

Alle bisher diskutierten Ergebnisse sind ein Indikator für den zusätzlichen Informations- 
gehalt, der bei einer konsequenten Berücksichtigung der Rauigkeit beider Kontaktkör- 
per aus den vorliegenden Simulationen gewonnen werden kann. Grundsätzlich lassen 
sich unabhängig vom Betriebspunkt für alle isotropen Konfigurationen in Tabelle 4.1 
qualitativ gleiche Ergebnisse finden. Die Resultate aller 231 simulierten Betriebspunkte 
werden daher für die drei Konfigurationen Isotrop 1-3 kompakt in den Abbildungen 4.5 
und 4.6 zusammengefasst. 

Die Erwartungswerte für die Scherspannungen E [TA] sind in den Abbildungen 4.5(a), 
4.6(a) und (e) als Funktionen der Erwartungswerte für den Druck E [p 4] dargestellt. Jede 
Kurve stellt die Erwartungswerte bei einer konstanten Gleitgeschwindigkeit vre und 
verschiedenen Einrückungen dar, wie sie zum Beispiel in Abbildung 4.4(b) zu sehen 
sind. Die eingezeichneten Orientierungslinien für den (deterministischen) Reibwert 
Ya sind dabei nur in guter Näherung mit E[uA] gleichzusetzen. Die Erwartungswerte 
für den Reibwert E[uA] sind in Tabelle 4.3 aufgelistet. Die Erwartungswerte für die 
Temperaturen E [04] sind in gleicher Weise in den Abbildungen 4.5(b), (b) und (f) als 
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Abbildung 4.4: Ausgewählte Ergebnisse für die eingelaufene Konfiguration Isotrop 1 bei einer Gleitgeschwin- 
digkeit vre] = 0.5 m/s und einer Einrtickung d = 2.08 um 
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(c) Verschiedene Autokovarianzfunktionen für die (d) Verteilungsdichtefunktionen pa, ,(Ay) der 
neuen Oberflachen (links) und die eingelaufenen Reibwertschwankungen 
Oberflachen (rechts) 


Abbildung 4.5: Zusammenfassende Ergebnisse fiir die Konfiguration Isotrop 1 


Funktionen der Erwartungswerte für den Druck E [p4] gezeigt. In allen Fällen markieren 
schwarze Linien den Fall Rg = 0 MPa/K und graue Linien den Fall Rọ = —1 MPa/K. 
Des Weiteren markieren gestrichelte Linien die Ergebnisse fiir die neuen Oberflachen 
im rein elastischen Lastfall, für den die Harte H = œ gewählt wird. 

Im rein elastischen Belastungsfall nimmt die Scherspannung E [TA] proportional zum 
Druck E [p4] zu, obwohl die Werte von Schiefe und Kurtosis von der äquivalenten Ober- 
fläche variieren (vgl. Abbildungen 4.1(b), (d) und (f)). Durch die Proportionalitat lässt 
sich dadurch auch für den gleitenden Kontakt zweier Körper mit rauen Oberflächen bei 
einer vergleichbaren tatsächlichen Kontaktfläche ein lastunabhängiger Reibwert folgern. 
Der niedrige Reibwert kann auf eine kleine tatsächliche Kontaktfläche bei vergleichs- 
weise hohen Drücken zurückgeführt werden. Die thermoelastischen Verschiebungen 
scheinen vernachlässigbar zu sein, während Rg den Reibwert geringfügig verringert. 
Durch die Proportionalität zwischen tatsächlicher Kontaktfläche und Normallast steigt 
auch die Temperatur E [04] linear mit dem Druck E [p4] an. 

Im Vergleich zum rein elastischen Lastfall hängt die Scherspannung E [t4] bei den 
eingelaufenen Oberflächen nichtlinear mit dem Druck E [p4] zusammen. Dies ist die 
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(h) Verteilungsdichtefunktionen fiir Isotrop 3 


Abbildung 4.6: Zusammenfassende Ergebnisse fiir die Konfigurationen Isotrop 2 und Isotrop 3 
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Folge der linksschiefen Oberflächenhöhen. In Untersuchungen zu statischen, nicht- 
normalverteilten äquivalenten Oberflächen lassen sich qualitativ gleiche Resultate 
finden [52, 125, 131]. Die tatsächliche Kontaktfläche wächst zunächst vergleichsweise 
stark mit dem Druck E[pı] an, dies führt zu einem hohen Reibwert bei niedriger 
Belastung. Der Einfluss der Gleitgeschwindigkeit und folglich der thermoelastischen 
Verschiebungen auf den Reibwert fällt größer aus, die Berücksichtigung einer tempe- 
raturabhängigen Versagensscherspannung durch Rg verringert den Reibwert deutlich 
stärker. Die Temperatur E [04] steigt aufgrund der nichtlinear wachsenden tatsächlichen 
Kontaktfläche ebenfalls nichtlinear mit dem Druck E [pa] an. 

Außerdem sind die gemittelten Autokovarianzfunktionen des Drucks 


J 
1 (i) 
Kpapa(Ax) = T SS Kýapa (Ax), 
j=l 
der Scherspannung 
Lé vi 
KraralOx) = 7 A. Kiara (An) 
j=l 
und des Reibwerts 
14 vii 
Eu, (Ax) = 7 ES Ee (Ax) 
j=l 


in den Abbildungen 4.5(c), 4.6(c) und (g) gezeigt, wobei hier die Ergebnisse von allen 
simulierten Betriebspunkten genutzt werden. Um einen direkten Vergleich zwischen 
den neuen und den eingelaufenen Oberflächen zu ermöglichen, beinhalten die Abbil- 
dungen die Verläufe für beide Lastfälle. Neben der durch die periodischen Oberflächen 
verursachte Kreisfrequenz o, = nt mm! sind offensichtlich auch weitere Frequenzan- 
teile enthalten und die Verläufe der Kurven unterscheiden sich in Abhängigkeit der 
Konfiguration sichtbar. Die ablesbaren Korrelationslängen Ar: für die Reibkraft liegen 
in derselben Größenordnung wie die von Rabinowicz aus Experimenten ermittelte 
Korrelationslänge AXxorr X 0.013 mm [215]. Die jeweiligen Korrelationslängen AXkorr 
für den Reibwert sind in Tabelle 4.3 aufgelistet. 

Die Abbildungen 4.5(d), 4.6(d) und (h) beinhalten die Verteilungsdichtefunktionen 
Dau, (Au) der Reibwertschwankungen 


Aua = pa -E [ua], 
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wobei hierfür ebenfalls die Ergebnisse von allen simulierten Betriebspunkten genutzt 
werden. Durch diese Vorgehensweise gehen etwaige Abhängigkeiten der Verteilungs- 
dichtefunktion Dau, (AU) von der Last verloren, als erster Ansatz soll dies jedoch genü- 
gen. Auch hier markieren die schwarze Linien den Fall Re = 0 MPa/K, graue Linien den 
Fall Ro = -1 MPa/K und gestrichelte Linien symbolisieren den rein elastischen Lastfall. 
Es ist erkennbar, dass die Standardabweichung für den rein elastischen Fall jeweils 
kleiner ausfällt. Die Standardabweichungen (E[A 4D" ? für die Reibwertschwankungen 
sind in Tabelle 4.3 aufgelistet. 


4.2  Anisotrope Oberflächen 


Die Auswertungen werden für die anisotropen Konfigurationen aus Tabelle 4.1 wieder- 
holt, um Unterschiede sowie Gemeinsamkeiten herauszuarbeiten. 

Die Verteilungsdichtefunktionen pz(z) der Oberflächenhöhen sind in den Abbildun- 
gen 4.7(a), (c) und (e) für den neuen und eingelaufenen Zustand der anisotropen 
Konfigurationen aus Tabelle 4.1 dargestellt. In allen drei Fällen sind die Oberflächenhö- 
hen nach dem Einlaufprozess zwar ebenfalls linksschief Ssx < 0 und die Werte für die 
Standardabweichung S, sowie die Kurtosis Sku haben sich verringert, allerdings fallen 
die Veränderungen im Vergleich zu den isotropen Konfigurationen deutlich kleiner aus. 
Die erhaltenen Ergebnisse passen dennoch auch hier gut zu Experimenten bezüglich 
dem Einlaufverhalten rauer Körper mit anisotropen Oberflächenstrukturen [124, 233], 
obwohl die aktuelle Modellierung ein vergleichsweise einfaches Verschleißverhalten 
beinhaltet. Ein Vergleich des Verschleißvolumens offenbart zudem, dass das Verschleiß- 
volumen von Konfiguration Anisotrop 1 über Anisotrop 2 zu Anisotrop 3 abnimmt, 
während es bei den isotropen Konfigurationen am größten ausfällt. Dies deckt sich mit 
Untersuchungen zum Verschleißvolumen bei isotropen, transversalen und lateralen 
Oberflächenstrukturen [233]. Grundsätzlich kann der geringere Verschleiß auf die 
strukturierteren Oberflächenrauigkeiten der anisotropen Konfigurationen zurückge- 
führt werden, die in diesem Fall lokal kleinere Lastspitzen verursachen. 

Des Weiteren sind die Werte von Kurtosis Sku und Schiefe S,; in den Abbildungen 4.7(b), 
(d) und (f) für die jeweiligen äquivalenten Oberflächen aus Gleichung (3.62) während 
einer Gleitperiode von Ax = 2 mm sowohl für die neuen als auch für die eingelaufenen 
Oberflächen dargestellt. Analog zu den isotropen Konfigurationen zeigen die Verläufe 
auch hier, dass die Oberflächenhöhen der äquivalenten Oberfläche selbst bei normalver- 
teilten Oberflächenhöhen der Kontaktkörper nicht als normalverteilt modelliert werden 
sollten. Zudem streuen die Werte insbesondere bei Konfiguration Anisotrop 1 deutlich 
stärker. Die Erwartungswerte von Standardabweichung E[S,], Kurtosis E[Sx„] und 
Schiefe E[S;x] der jeweiligen äquivalenten Oberflächen sind für die Gleitperiode von 
Ax = 2 mm ebenfalls in Tabelle 4.3 aufgelistet. 
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Abbildung 4.7: Ausgewählte Eigenschaften der neuen und eingelaufenen Oberflächenhöhen von den 
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Konfigurationen Anisotrop 1, Anisotrop 2 und Anisotrop 3 aus Tabelle 4.1 
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Abbildung 4.8: Bereiche auf den Oberflächen der drei anisotropen Konfigurationen, die mindestens in einem 


inkrementellen Zeitschritt in Kontakt stehen bei der Gleitgeschwindigkeit oi = 0.5 m/s und 
vergleichbarer Last 
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Erneut werden die neuen und eingelaufenen Oberflächen für die folgenden Simulati- 
onsergebnisse verwendet und die Simulationen bis zum Erreichen eines stationären 
Betriebszustandes durchgeführt. Die erhaltenen und nachfolgenden Ergebnisse sind 
daher periodisch. Für die einzelnen Betriebspunkte können Blitztemperaturen, thermo- 
elastische Verformungen und oszillierende kontaktmechanische Größen aufgrund der 
zweiseitigen Rauigkeit gefunden werden. 

Die Ergebnisse spiegeln dabei die anisotropen Oberflächenstrukturen wider. Dies wird 
anhand von den Abbildungen 4.8(a)-(f) deutlich, in denen für alle anisotropen Konfi- 
gurationen die Bereiche der Oberflächen in schwarz markiert sind, die während des 
gesamten Gleitwegs von Ax = 2 mm mindestens einmal in Kontakt stehen. Wiederrum 
sind deutlich mehr Oberflächenbereiche am eigentliche Reibprozess beteiligt, als rein 
statische Simulationen vermuten lassen würden. 

Die kompakte Zusammenfassung der Ergebnisse aller simulierten Betriebspunkte 
erfolgt in den Abbildungen 4.9 und 4.10. Im Vergleich zu den isotropen Konfigura- 
tionen fällt in den Abbildungen 4.9(a), 4.10(a) und (e) der größere Reibwert bei rein 
elastischen Belastungsfällen auf, der auf eine zwar immer noch proportional, aber nun 
schneller mit dem Druck E [p4] ansteigende tatsächliche Kontaktflache zurückgeführt 
werden kann. Dies resultiert in den Abbildungen 4.9(b), 4.10(b) und (f) in ebenfalls 
höheren Temperaturen E [O4]. Ansonsten ähneln die Verläufe der Erwartungswerte 
von Druck E [pa], Scherspannung E[t,] und Temperatur E [04] weitestgehend den 
Verläufen für die isotropen Konfigurationen, sodass auch die Schlussfolgerungen gleich 
bleiben. Unterschiede im Detail lassen sich auf die verschiedenen Oberflächenpaarun- 
gen zurückführen. 

Anders sieht es bei den Autokovarianzfunktionen in den Abbildungen 4.9(c), 4.10(c) 
und (g) aus. Hier sind die Auswirkungen der anisotropen Leistungsdichtespektren 
zu beobachten, die vor allem bei Konfiguration Anisotrop 3 in Abbildung 4.10(g) zu 
einer großen Korrelationslänge führen. Im Vergleich zu den isotropen Konfigurationen 
besitzen die erhaltenen Verteilungsdichten in den Abbildungen 4.9(d), 4.10(d) und (h) 
außerdem eine höhere Varianz. Durch das geringere Verschleißvolumen, ergeben sich 
dort kleinere Unterschiede zwischen den neuen und eingelaufenen Oberflächenpaarun- 
gen. Die jeweiligen Korrelationslängen AxXxorr und Standardabweichungen (lau Di 2 
der berechneten Reibwerte werden in Tabelle 4.3 aufgelistet. 
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Abbildung 4.9: Zusammenfassende Ergebnisse für die Konfiguration Anisotrop 1 


(d) Verteilungsdichtefunktionen pa, (Au) der 
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Abbildung 4.10: Zusammenfassende Ergebnisse für die Konfigurationen Anisotrop 2 und Anisotrop 3 
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4.3 Einordnung der Ergebnisse 


Die durchgeführten Simulationen ermöglichen bisher neuartige Einblicke in den Kon- 
takt gleitender Körper. Die Berücksichtigung der Rauigkeit beider Kontaktkörper 
offenbart dabei transiente Eigenschaften, die durch Modellierungen mit einer einsei- 
tigen Rauigkeit nicht erfasst werden können. Hierzu zählen in besonderem Maße 
die Blitztemperaturen und die transiente Kontaktfläche, die wiederum transiente 
Spannungsfelder implizieren. Es zeigt sich ein starker Einfluss des Einlaufprozesses auf 
die Ergebnisse, wobei der Verschleiß bei den anisotropen Konfigurationen vergleichs- 
weise geringer ausfällt. Durch den Einlaufprozess erhöht sich außerdem bei allen 
Konfigurationen der deterministische Anteil des Reibwerts. Folglich können bereits 
hier reibungserregte Schwingungen entstehen, die wiederum durch das Minimalmodell 
der nicht-konservativen Kopplung aus Abschnitt 2.3.2 erklärbar sind. 

Zusammenfassend lässt sich für alle sechs Oberflächenkonfigurationen aus Tabelle 4.1 
sowohl eine Last- als auch eine Geschwindigkeitsabhängigkeit des deterministischen 
Reibwertanteils feststellen. Die Lastabhängigkeit des Reibwerts tritt nur bei den Ober- 
flächen mit E[S,,] + 0 auf, während sich für die normalverteilten Oberflächen mit 
ESA ~ 0 auch bei einer Relativbewegung der Kontaktkörper ein lastunabhängiger 
Reibwert ergibt. Letzterer ist auf eine Proportionalität zwischen der tatsächlichen 
Kontaktfläche und der Last aufgrund einer äquivalenten Oberfläche mit im Mittel 
immer noch normalverteilten Höhenwerten zurückzuführen (vgl. Tabelle 4.3). Ob 
eine Lastabhängigkeit vorliegt, lässt sich somit anhand der Schiefe E[S;x] beurteilen. 
Die Geschwindigkeitsabhängigkeit des Reibwerts ist einerseits durch die Abhängig- 
keit der Versagensscherspannung von der Gleitgeschwindigkeit und der Temperatur 
sowie andererseits durch die Berücksichtigung der thermoelastischen Ausdehnungen 
begründbar. Insbesondere die bereits bei kleinen Relativgeschwindigkeiten auftre- 
tenden hohen Blitztemperaturen sind dabei für die Absenkung des Reibwerts von 
entscheidender Bedeutung. Grundsätzlich passen alle erhaltenen Ergebnisse zumindest 
qualitativ gut zu Ergebnissen in der Literatur, wobei sich die Simulationen durch die 
hergeleiteten Einflussfunktion leicht auf beliebige Materialparameterkombinationen 
erweitern lassen. Die berechneten Autokovarianz- und Verteilungsdichtefunktionen 
widersprechen einer Modellierung des stochastischen Reibwertanteils als normalver- 
teiltes weißes Rauschen. Diese Schlussfolgerung passt gut zu veröffentlichten Mes- 
sungen [215, 229]. Allerdings lässt sich in der Literatur noch ein generelles Defizit 
hinsichtlich gezielter Auswertungen des stochastischen Reibwertanteils festhalten. Die 
anisotropen Eigenschaften der Oberflächen sind in den Autokovarianzfunktionen des 
Reibwerts wiederzuerkennen. Darüber hinaus werden keine Einflüsse von Temperatur 
oder thermoelastischen Verformungen auf den stochastischen Anteil des Reibwerts 
gefunden. Ob die verwendete Diskretisierung für allgemeingültige Aussagen ausreicht, 
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Tabelle 4.3: Zusammenfassung der Ergebnisse aller isotropen und anisotropen Konfigurationen (gerundet 
auf die dritte Nachkommastelle) 


Konfiguration Isotrop 1 Isotrop 2 Isotrop 3 
Harte co 3Re o0 3Re co 3Re 
E[S,] in um 0.987 | 0.853 | 0.991 | 0.835 0.99 0.806 
E[Sku] 2.91 2.839 2.99 2.905 2.98 2.973 
E[Ssx] —0.005 | -0.376 | 0.004 | -0.455 | 0.003 | -0.527 


Amen In mm | 0.077 | 0.08 | 0.076 | 0.103 | 0.048 | 0.069 
Tog /4 in mm 1/16 | ~1/16| 1/24 | ~1/24| 1/32 | ~1/32 


Elup] 0.055 -= 0.041 u 0.034 = 
E[ua] 0.109 - 0.105 = 0.076 = 
(lau Di 0.008 | 0.038 | 0.008 0.03 0.004 | 0.025 
Konfiguration Anisotrop 1 Anisotrop 2 Anisotrop 3 
Harte co 3Re co 3Re oo 3Re 
E[S,] in um 0.979 | 0.887 | 0.993 | 0.907 | 0.992 | 0.931 
E[Sku] 2.913 | 2.838 | 3.053 | 2.879 | 2.804 | 2.713 
EIS —0.027 | —0.279 | —0.021 | —0.297 | —0.055 | —0.236 


AXkorinmm | 0.139 | 0.25 | 0.134 | 0.237 | 0.478 | 0.418 
T.«/4 in mm 1/16 | =1/16 | 3/32 | ~3/32 | 3/16 | ~3/16 
E[up] 0.071 u 0.071 = 0.071 — 
E[ua] 0.191 — 0.155 — 0.194 — 
(E[Au2 Di 0.026 | 0.072 | 0.012 | 0.057 | 0.029 | 0.143 


kann auch aufgrund einer fehlenden Validierung mit Experimenten nicht abschließend 
beantwortet werden. Almqvist [6] folgert jedoch aus seinen Untersuchungen, dass 
bereits bei einer Berücksichtigung von 25 % der tatsächlich vorhandenen Frequenzan- 
teile repräsentative Ergebnisse zu erwarten sind. Daher scheinen die Aussagen des 
Modells im Rahmen der Halbraumannahme gültig zu sein. 

Es verbleibt die Frage, inwieweit sich die erhaltenen Ergebnisse auf die Oberflächenrau- 
igkeiten zurückführen lassen. Bezüglich des deterministischen Reibwerts wird diese 
Fragestellung zumindest für den quasi-statischen Normalkontakt bei normalverteil- 
ten Oberflächenhöhen bereits weitestgehend durch vorhandene Veröffentlichungen 
behandelt und lässt sich durch die Vernachlässigung von thermischen Effekten auf 
gleitende Kontakte erweitern. Bush, Gibson und Thomas [41] sowie Persson [201] folgern 
für vergleichsweise kleine tatsächliche Kontaktflächen eine Proportionalität zwischen 
der tatsächlichen Kontaktfläche und dem Kehrwert des quadratischen Mittelwertes 
der Oberflächenhöhengradienten (E[||Vz||2])~!/2. Daraus lässt sich ebenfalls auf eine 
Proportionalität zwischen (E[]|Vz]]?])-'/? und dem Reibwert u4 beziehungsweise der 
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Temperatur 04 schließen, da beide proportional zur tatsächlichen Kontaktflache sind. 
Mit dem Proportionalitätsfaktor ¥87 von Persson [201] kann der Reibwert 


1-w 1-7? 
bn = Vir 14 | Ts 


mE, = Ee | VENIVA 


hergeleitet werden, wobei die zugehörigen Erwartungswerte E[up] über eine Gleit- 
periode von Ax = 2 mm in Tabelle 4.3 aufgelistet sind. Um dabei die numerisch 
fehleranfällige Berechnung der Gradienten über finite Differenzen zu vermeiden, 
werden stattdessen die spektralen Momente der jeweiligen äquivalenten Oberflächen 


1 
ELIF] = 25 |] S2z@lol?do 
TI 
R2 


ausgewertet [178]. Der Korrelationskoeffizient p(E[ur], EluA]) ~ 0.929 fällt zwar erwar- 
tungsgemäß hoch aus, sollte aber wegen der geringen Anzahl an Stichproben noch durch 
weitere Simulationen bestätigt werden. Der absolute Unterschied in den Reibwerten 
E[up] und E[uA] lässt sich zudem auf mehrere Faktoren zurückführen. Zum einen 
beziehen sich die Herleitungen von Persson [201] aufnormalverteilte Oberflächenhöhen, 
wohingegen die Auswertungen in den Abbildungen 4.1 und 4.7 zeigen, dass dies bei 
gleitenden Konfigurationen nicht durchweg der Fall ist. Zum anderen unterschätzt 
Persson die tatsächliche Kontaktfläche [188], während die für die Simulationen gewählte 
Diskretisierung sicherlich zu einer Überschätzung der tatsächlichen Kontaktfläche führt. 
Für die eingelaufenen Oberflächen entfällt die Proportionalität vollends, sodass für 
einen Zusammenhang zwischen den charakteristischen Oberflächenwerten und dem 
Reibwert auf numerische Auswertungen zugegriffen werden muss. 

Für einen Zusammenhang der erhaltenen Korrelationslängen des Reibwerts zu den 
Oberflächenstrukturen wird zunächst die äquivalente Oberfläche zweier Körper mit 
den harmonischen Oberflächenhöhen zı = Asin(wıx) und z2 = Bsin(w2x) durch 
Gleichung (3.62) 


z(x,t) =Z1(x) + z2(x — Vret ) 


+ - 
=(A + B)sin (= 7 o? y ell cos (= S ah Ze) 
w1 + w2 W20 re] ` (= - (an W20 re] 
= z t 
+(A B) cos ( SH > t > * 5 


berechnet. Hieraus können direkt die effektiven Periodendauern Teg = 47/|w 1 + w2| 
abgelesen werden. Mit der kleineren der beiden Periodendauern folgt die dominierende 
Korrelationslange Axa = Terr/4, wobei die Amplituden keinen Einfluss zeigen. Die auf 
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diese Weise hergeleitete Korrelationslange wird in Tabelle 4.3 fiir alle Konfigurationen 
aus Tabelle 4.1 aufgelistet, wobei on und ous mit den jeweiligen Werten von wr 
identifiziert werden. Bei den anisotropen Konfigurationen ist dabei jeweils noch 
Yy und Wy zu berücksichtigen. Der Korrelationskoeffizient P(AXkorr, Terr/4) = 0.934 
fällt ebenfalls hoch aus, sodass weiterführende analytische Untersuchungen hierauf 
aufbauen können. 

Eine vergleichbar hohe Korrelation konnte für die Reibwertschwankungen E[A A nicht 
gefunden werden, weshalb eine endgültige Abschätzung der Schwankungsbreite aus 
den beteiligten Oberflächen offen bleibt. 


Da die bisherigen Simulationsergebnisse lediglich eine Temperaturabhängigkeit der Ver- 
sagensscherspannung und der Härte berücksichtigen, wird abschließend eine Erweite- 
rung auf möglichst allgemeine temperaturabhängige Materialparameter diskutiert. Die 
zugrundeliegenden Gleichungen für temperaturabhängige Wärmeleitfähigkeit K(T), 
spezifische Wärmekapazität ca(T) und Wärmeausdehnungskoeffizient a(T) lauten [181] 


T 
DEE EE 
To 


d 
V- (K (T) VO) = pea (T) Ê. 


Durch die Transformation [181] 


T 


d K(9) 
T=To+ d Grat (4.1) 


To 


können diese Gleichungen zu 


uAr + (A + u) V (V - r) = DÄ + 2u) ott, (4.2) 
or 
K(T)AT = = (4.3) 


umformuliert werden, wobei von einer Proportionalitat zwischen a(T) und K(T) ausge- 
gangen wird. Es folgt weiterhin die Beziehung zwischen Temperatur und Warmestrom 


q = —K(T)V@ = —K(To)VI, 


während x(T) = K(T)/pca(T) die temperaturabhängige Temperaturleitfahigkeit bezeich- 
net. Die Gleichungen (4.2) und (4.3) stimmen strukturell mit den homogenen Gleichun- 
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gen (2.15) und (2.16) für den einseitig gekoppelten Fall mit € — 0 überein. Da für viele 
Materialien in guter Näherung von einer temperaturunabhängigen Temperaturleitfä- 
higkeit x(T) ~ konst. ausgegangen werden kann [15], ist die Verwendung der bereits 
in Unterkapitel 3.1 hergeleiteten Einflussfunktionen mit nun neuer Interpretation, 
zumindest im Rahmen der getroffenen Annahmen, weiterhin ohne Probleme möglich. 
Allerdings müssen dafür die Parametrisierungen von K(T), ca(T) und a(T) explizit 
bekannt sein. 

Im Folgenden werden die Auswirkungen von temperaturabhängigen Materialparame- 
tern in der Form 


K(T) ~K(To)(1 + y(T - To)), 
a(T) ~a(To)(1 + y(T - To)) 


untersucht, wobei Wärmeleitfähigkeit K(T) und Warmeausdehnungskoeffizient a(T) 
jeweils durch eine Taylorreihe erster Ordnung approximiert sind. Eine solche Parame- 
trisierung bietet sich im Zusammenhang mit der zugrundeliegenden Annahme für 
die Wärmeleitungsgleichung |0/To| « 1 an. Außerdem wird entsprechend von einer 
konstanten Temperaturleitfähigkeit x(T) = K(T)/pca(T) ~ konst. ausgegangen. Die 
Gleichung (4.1) lässt sich zu 


1 
T-h= ant -hy +T-h 


auswerten, es verbleibt nur eine physikalische Lösung 


2(T — To) 
14/1 +270 = 1) 


Für die thermoelastischen Verschiebungen ergeben sich wegen der Proportionalität von 
Wärmeleitfähigkeit K(T) und Warmeausdehnungskoeffizient a(T) keine Unterschiede. 

Anhand eines parabolischen Wärmestroms q = qo (1 — r?/ r2)" * innerhalb des Gebietes 
r<fqmitgqo =1 MJ/sm?, ra = 0.5 mm und K(To) = 40 J/msK werden die Auswirkun- 
gen verdeutlicht. In den Abbildungen 4.11(a) und (b) ist die relative Temperaturdifferenz 
0 - O|,-9/ O|,-9 in % für die beiden Werte y = +0.01 KI entlang der r-Richtung 
und der dimensionslosen Zeit t = «t/r? dargestellt. Die stationären Werte sind als 


0=T-N= 


gestrichelte Linien ebenfalls eingezeichnet. Der Einfluss der temperaturabhängigen 
Materialparameter ist deutlich zu erkennen und es lässt sich feststellen, dass die 
Berücksichtigung der Temperaturabhängigkeit durchaus eine Rolle spielt. Da bei 
vielen Materialien y > O gilt, mindert die Berücksichtigung von temperaturabhängigen 
Materialparameter zudem die Auswirkungen einer mit der Temperatur sinkenden 
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-6 4 > 
0 25 50 75 100 
r/Ya xt /12 
(a) Relative Temperaturdifferenz 0 — O|,=0 j d (b) Zeitlicher Verlauf der relativen Temperaturdiffe- 
in% firy = +0.01 Kl zu verschiedenen renz 0 — du / HN in % für y = £0.01 KI in 
Zeitpunkten entlang der r-Richtung der Mitte der Wärmequelle 


Abbildung 4.11: Einfluss von temperaturabhangigen Materialparametern auf das Temperaturfeld 


Versagensscherspannung. Die zusätzliche Konstante y bietet darüber hinaus eine 
Möglichkeit, die Simulationen an Experimente anzupassen. Eine Reihe verschiedener 
weiterer Parametrisierungen können der Literatur entnommen werden [166, 181]. Da 
die Ergebnisse der thermomechanischen Kontaktsimulationen auf einen eher gerin- 
gen Einfluss der thermoelastischen Verschiebungen auf kontaktmechanische Größen 
hindeuten, wird hier auf weiterführende Untersuchungen mit nicht-proportionalen 
Wärmeleitfähigkeiten K(T) und Wärmeausdehnungskoeffizienten a(T) verzichtet. 
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5 Diskussion im Kontext von 
reibungserregten Schwingungen 


Die thermomechanischen Kontaktsimulationen in Kapitel 4 haben insbesondere vier 
wesentliche Eigenschaften des berechneten Reibwerts aufgezeigt. Neben der Abhängig- 
keit des deterministischen Reibwertanteils von der Gleitgeschwindigkeit und der Last 
wurden auch genauere Erkenntnisse zur Verteilungsdichte und zur Autokorrelations- 
funktion des stochastischen Reibwertanteils erhalten. Das vorliegende Kapitel wird die 
Auswirkungen dieser Eigenschaften im Kontext von reibungserregten Schwingungen 
diskutieren. Dabei sind die Folgen der Geschwindigkeitsabhängigkeit bereits ausführ- 
lich in der Literatur behandelt (vgl. Abschnitt 2.3.1), jedoch fehlen Untersuchungen 
bezüglich Lastabhängigkeit und stochastischem Reibwert. Diese erfolgen jeweils in den 
Unterkapiteln 5.1 und 5.2. Da die Kontaktsimulationen zudem eine Modellierung mit 
zwei Freiheitsgraden nahelegen, wird dafür das Minimalmodell aus Abschnitt 2.3.2 
verwendet. Interessanterweise weist Ibrahim [110] bereits auf das Fehlen von For- 
schungsbemühungen hinsichtlich stochastischer Modellierungen bei Minimalmodellen 
mit nicht-konservativer Kopplung hin, entsprechende Arbeiten konnten aber nicht 
gefunden werden. 


5.1 Lastabhängiger Reibwert 


Für einen last- und damit einrückungsabhängigen deterministischen Reibwert ua(z) 
verändern sich die gekoppelten Differentialgleichungen für den Zustand des Gleitens 
aus Gleichung (2.39) in Abschnitt 2.3.2 zu 


m O)\ [x A dı O\ {x a c11 C12 x\ _ (UalZ)enz sign(drei) 
0 m]\ä 0 do} \z Cou C2 + Cy} \Z u Fe j 
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die sich um die Ruhelage (xo, Zo)" linearisieren lassen 


m 0 Ax + dı 0 Ax 4 C11 C12 — Du Ax _ 0 
0 m] \Az 0 do} \Az co «Con +n | \Az} ol? 


Ou 
worin fi = ua(Zo) + zo Déi 
Z=Z 


bezeichnet. Für beliebige funktionale Zusammenhänge 


Hatz) lässt sich die Ruhelage (xo, zo)" in der Regel nur noch numerisch berechnen, die 
Existenz mehrerer Ruhelagen ist möglich. Im Vergleich zu den linearisierten Gleichun- 
gen (2.40) in Abschnitt 2.3.2 ergibt sich ein zusätzlicher Einfluss des Reibwertgradienten 
auf die Stabilität der Ruhelage, sodass auf eine oszillatorische Instabilität der Ruhelage 
in Abhängigkeit von fi geschlossen werden kann. Das qualitative Stabilitätsverhalten 
bleibt allerdings aufgrund der unveränderten Struktur der Gleichungen erhalten. 


5.2 Stochastischer Reibwert 


Für einen stochastischen Reibwert M, wird eine möglichst allgemeine Parametrisierung 
von Verteilungsdichte- und Autokorrelationsfunktion gewählt, sodass anschließende 
Parameterstudien erfolgen können. Um die Untersuchungen auf die Auswirkungen 
des stochastischen Reibwertanteils zu beschränken, wird im Folgenden von einem 
konstanten deterministischen Reibwertanteil ua = EIM;] = konst. ausgegangen. Die 
Diskussion erfolgt in zwei Schritten. Zunächst wird ein auf Ortsebene beschriebener 
stochastischer Prozess in Abschnitt 5.2.1 derart modelliert, dass er einerseits die Resul- 
tate der thermomechanischen Simulationen abbildet und andererseits physikalisch 
konsistent ist. Anschließend erfolgt die Transformation in den Zeitbereich, um eine 
Verträglichkeit mit den Differentialgleichungen des nicht-konservativen Minimalm- 
odells herzustellen. Die Auswirkungen des auf diese Weise modellierten stochastischen 
Reibwerts auf die Dynamik des Systems erfolgt daraufhin in Abschnitt 5.2.2. 


5.2.1 Reibwertmodellierung 


Für die Autokorrelationsfunktion wird im Folgenden die Parametrisierung 


Rusu _ x 
Ge = e7POxlAx| (cos EI VI=D2jaxl) + = —— sin Ei V1- GË (5.1) 


gewählt, was dem zweiseitigen Leistungsdichtespektrum 


Bib, _ 4DO3 


K (Q2 - w2)? +4D20202 
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entspricht. Beide Ausdrücke beinhalten die Filterparameter Q, und D, die für die 
Anpassung an experimentell oder rechnerisch erhaltene Ergebnisse herangezogen 
werden können. Darüber hinaus bezeichnet R, hier und im Folgenden die Standardab- 
weichung des stochastischen Reibwerts. Für die Verteilungsdichte soll die Modellierung 
physikalisch konsistent und daher beschränkt erfolgen, sodass ein normalverteilter 
Prozess ausgeschlossen ist. 

Zhu und Cai [266] leiten die stochastischen Filtergleichungen im Sinne Itös 


0 1 /4D03 0 
Xx = X,dx + 4/ ———= 
SS Si ES MOON Dee ee 


her, die die frei wählbaren Filterparameter Qx, D, A und ı beinhalten. Wird der 
stochastische Reibwertanteil mit us = Xı identifiziert, ergibt sich für den stationären 


dw, (5.2) 


Fall die symmetrische und auf das Intervall X; € [-A, A] beschränkte Verteilungsdich- 
tefunktion 


ERR 
X, 
A2 


T(2ı +2) 
2241 to +1)? 


Pus(X1) = H (5.3) 


die lediglich von der Gamma-Funktion I...) und den beiden Filterparametern A 
sowie ı abhängt. Es folgen der Mittelwert E[M,] = ua, die Varianz E IM; - La)? | = 
A? /(21 + 3) und beliebige höhere zentrale Momente 


A” (1+ (-1)")T 


E[(M; = ua)" ] = (5.4) 


Des Weiteren besitzen die stationären Lösungen Xı von Gleichung (5.2) allesamt die 
Autokorrelationsfunktion (5.1). 

Die gewählte Autokorrelationsfunktion aus Gleichung (5.1) wird in Abbildung 5.1(a) mit 
der in Unterkapitel 4.1 berechneten Autokovarianzfunktion von Konfiguration Isotrop 2 
aus Abbildung 4.6(c) verglichen. Der Übergang von den periodischen Ergebnissen auf 
nicht-periodische Oberflächenkonfigurationen erfolgt hierbei durch die Simulation mit 
den Filtergleichungen automatisch. Darüber hinaus wird in Abbildung 5.1(b) auch die 
Verteilungsdichtefunktion aus Gleichung (5.3) mit der in Unterkapitel 4.1 berechneten 
Verteilungsdichtefunktion von Konfiguration Isotrop 2 aus Abbildung 4.6(d) verglichen. 
Beide Parametrisierungen scheinen die erhaltenen Verläufe aus den thermomechani- 
schen Kontaktsimulationen zufriedenstellend abzubilden. Durch die Berücksichtigung 
weiterer Frequenzanteile wird eine noch bessere Übereinstimmung erwartet, aus 
Gründen der Übersichtlichkeit erfolgt hier jedoch die Beschränkung auf die einzelne 
Kreisfrequenz Qx. 
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70 
— Pu; (Ap) 


35 


-1 -0.5 0 0.5 1 -0.1  -0.05 0 0.05 0.1 
Ax inmm Au 
(a) Autokovarianzfunktion aus Gleichung (5.1)im (b) Verteilungsdichtefunktion aus Gleichung (5.3) im 
Vergleich zu Abbildung 4.6(c) Vergleich zu Abbildung 4.6(d) 


Abbildung 5.1: Vergleich von Autokorrelations- und Verteilungsdichtefunktion mit den erhaltenen Ergeb- 
nissen der thermomechanischen Kontaktsimulation für Konfiguration Isotrop 2 in Unterka- 
pitel 4.1 


Durch den Zusammenhang dx = |v,.ıldt können die Filtergleichungen vom Orts- in 
den Zeitbereich transformiert werden. Um danach immer noch einen entkoppelten 
stationären Prozess zu erhalten, soll an dieser Stelle von einer konstanten Relativge- 
schwindigkeit drei = |vo — ž| = vo ausgegangen werden. Die Umwandlung aus dem 
Orts- in den Zeitbereich 


dx i oea : dw, 55 
+ S 
"Lo ap0ol"" 2ı+1\)2-xX2-x2/02] ° 


kann dann mit Q; = O;v9 problemlos erfolgen. Die Verteilungsdichtefunktion bleibt 
unverändert, dementsprechend sind Mittelwert und Standardabweichung unabhängig 
von der Geschwindigkeit. Allerdings lässt sich jeweils ein Einfluss der Geschwindig- 
keit vo auf die Autokorrelationsfunktion 


Riva D 
et: _ e DO [os (ovi = D3!) + —— sin (avi = rell 
Vi-D 


2 
R3 
und das zweiseitige Leistungsdichtespektrum 


Sonn, 4DQ? 
— = — (5.6) 
R3 (02 - w?) + 4D Qw? 


feststellen, sodass sowohl das Leistungsdichtespektrum für vo — 0 als auch die 
Korrelationslänge des Prozesses fiir ou — œ gegen null konvergiert. Die Auswirkungen 
lassen sich dementsprechend physikalisch plausibel interpretieren. Q; kann als systems- 
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0 4 
—1 —0.5 0 0.5 1 


u/A wi /Qi 
(a) Verteilungsdichtefunktion aus Gleichung (5.3) für (b) Zweiseitiges Leistungsdichtespektrum aus 
verschiedene Filterparameter ı Gleichung (5.6) für D = 0.1,0.3, v2/2 


Abbildung 5.2: Verteilungsdichtefunktion und Leistungsdichtespektrum des Reibwerts für verschiedene 
Filterparameter 


pezifische Kreisfrequenz aufgefasst werden, die im vorliegenden Fall durch die Oberflä- 
chenrauigkeiten verursacht wird. Sie ist auf natürliche Weise mit der Geschwindigkeit vo 
gekoppelt. Der Parameter D beeinflusst die Abklingrate der Autokorrelationsfunktion, 
während A und ı maßgebend für die Form der Verteilungsdichtefunktion sind. Sind 
die Werte für Standardabweichung R, und Kurtosis Rx, der zu approximierenden 
Verteilungsdichtefunktion aus Experimenten oder Simulationen bekannt, können die 
Filterparameter ı und A beispielsweise zu ı = (9 — 5Rx„)/(2Rxu — 6) und anschließend 
A = RyV2ı+3 gewählt werden. Diese Zusammenhänge lassen sich nach kurzer 
Rechnung aus Gleichung (5.4) herleiten. Dabei wird jedoch auch ersichtlich, dass der 
hier gewählte Prozess nur für Werte der Kurtosis im Bereich 3/2 < Rx. < 3 geeignet ist. 
Darüber hinaus lässt sich der Filterparameter Qy bei bekanntem Dämpfungsmaß D in 
Abhängigkeit der Korrelationslänge Ax or, durch Qy = 1/2AXkorrV1 — D? abschätzen. 
Ein direkter Zusammenhang zu den Oberflächenrauigkeiten der Kontaktkörper ist 
somit gegeben. 

Die Verteilungsdichtefunktion ist für verschiedene Werte von A und ı in Abbil- 
dung 5.2(a) dargestellt. Abbildung 5.2(b) zeigt den Verlauf des Leistungsdichtespek- 
trums für verschiedene Dämpfungsmaße D. Die Maxima der Leistungsdichtespektren 
bei w;/Q; = V1 -— 2D? sind durch die gestrichelte rote Linie verbunden. Die Abbil- 
dungen verdeutlichen, dass die vier Filterkoeffizienten Q;, D, A und ı vielseitige 
Anpassungsmöglichkeiten liefern und dadurch Parameterstudien ermöglichen. 

Um negative Werte unter der Wurzelfunktion der Filtergleichungen (5.5) während 
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den numerischen Simulationen zu vermeiden, wird unter Berücksichtigung von Itös 
Lemma [57, 227] auf die Transformation 


Xı =Asin(@) cos(0) 
X = - AQ; sin(¢) sin(0) 


zurückgegriffen [266]. Abschließend sei ergänzt, dass die Lösung der Filtergleichungen 
für die Anfangsbedingungen ¢ = n/2 und 0 = 0 sowie dem Dämpfungsmaß D = 0 in 
die harmonische Funktion Xı = Acos(Q;t) übergeht. 


5.2.2 Nicht-konservative Kopplung mit stochastischem 
Reibwert 


Für den stochastischen Reibwert M; = ua + us verändern sich die gekoppelten Differen- 
tialgleichungen für den Zustand des Gleitens aus Gleichung (2.39) in Abschnitt 2.3.2 


zu 
m O)\ (x e dı DI aon C12 x\ (ua + Us)enz sign(vo — X) 
0 m)\z 0 do} \z Cou C2 + Cy} \Z u Fe ý 


die sich um die Ruhelage (im Sinne der Erwartungswerte) 


UaCn—c12 
(2) — | cu(C22+cn)+c21(HaCn—c12) F 
= e 


Zo cu 
C11(C22+Cn )+C21 (Maen —C12) 


linearisieren lassen 

m 0O\ {Ax N dı 0\[Ax J Cu C12 — (Ha + Hs)Cn\ (Ax\ _ [HscnZo 

0 mila 0 d)/\A2) \ca C22 + Cn Az), E 
Hierbei wird von einer stochastischen Unabhängigkeit zwischen z und u, ausgegangen. 
Es ist erkennbar, dass der stochastische Reibwertanteil als Parameteranregung und 


Fremdanregung wirkt. Nach einer Umskalierung der Zeit 7 = t/to mit to = ym/cn 
ergeben sich die linearisierten Differentialgleichungen 


1 0\ (Ax” 2561 0 \ [Ax’ Cu Ce- (ua + us)\ (Ac) _ (sën 
( | en 0 A ee ee IC. 5 ) (5.7) 
me — S 


M D K+N+S(t) 


wobei die dimensionslosen Dampfungsbeiwerte ou = dr (Nie und 62 = da/2 mc, 
sowie die dimensionslosen Steifigkeitsparameter Cu = cı1/Ca, Ce = C12/¢n = C21/Cn 
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und Cy = C22/Cc, eingeführt werden. Zudem sind die Massenmatrix M = MI, die 
Dampfungsmatrix D = D', der symmetrische Anteil der Steifigkeitsmatrix K = KI, 
der zirkulatorische Anteil der Steifigkeitsmatrix N = -N" und der stochastische Anteil 
der Steifigkeitsmatrix S(t) mit E[S(z)] = 0 enthalten. Für die folgende Stabilitatsunter- 
suchung wird das homogene Gleichungssystem (5.7) betrachtet. 

Die Berechnung des größten Lyapunov-Exponenten 


Ina; 
Ay = lim : In a = lim l f dln As (5.8) 


für die triviale Lösung des homogenen Systems (5.7) lässt Stabilitätsaussagen mit Wahr- 
scheinlichkeit 1 zu. A; bezeichnet dabei die transiente Entwicklung der euklidischen 
Norm einer zum Zeitpunkt t = 0 vorliegenden Abweichung Ag von der Ruhelage. Wei- 
tere Informationen zu stochastischer Stabilitätstheorie können in Anhang D gefunden 
werden. Die Transformation in die neuen Koordinaten 


Ar =VAx2 + Ax? + Az? + Az”, 


tang, =Ax’/Ax, 
tan Y: =Az’/Az, 


tan 9, =VAz? + Az’? /V Ax? + Ax’? 


erlaubt unter Berücksichtigung von Itôs Lemma die Umformulierung des homogenen 
Gleichungssystems (5.7) in ein System von stochastischen Differentialgleichungen im 
Sinne von Itô 


din A; fi (Pr, Wr, Su, or, 9.) 0 
do, fr GE 0 
dy; = fa Pr, Wr, Sx, Pr, Kb 0 
ds, 7 fa (pr, Wr, EI or, 9.) deg 0 aWe, 69) 
der fs (Pr Or) fr (pr, Or) 
d0, Js (dr, Or) fs (Qr, d. 


wobei die ebenfalls in der Zeit neu skalierten Filtergleichungen (5.5) fiir den stochas- 
tischen Reibwert angehängt sind. Die vollständigen Ausdrücke für die Drift- und 
Diffusionsterme sind in Anhang D aufgelistet. Sie sind außerdem unabhängig von der 
Amplitude A.. Folglich hängt der größte Lyapunov-Exponent aus Gleichung (5.8) nur 
von den stationären Winkelprozessen Pr, Yr, 9x, Qr und 9, ab. 

Da das Lösen einer in diesem Fall fünfdimensionalen partiellen Differentialgleichung 
für die Verbundwahrscheinlichkeitsdichtefunktion immer noch mit großem Aufwand 
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verbunden ist [14, 247], wird auf eine numerische Zeitintegration (vgl. Abschnitt 2.1.3) 
mit zufälligen Anfangsbedingungen für die ersten vier Zustände In Ap ~ U(0,1), Po ~ 
U(0,2r), Po ~ U(0,27) sowie än ~ U(0,2r) zurückgegriffen, während für die 
Filtergleichungen du = n/2 und On = 0 verwendet wird. Die Simulation ist auf 
107 Zeitschritte mit einer Zeitschrittweite von At = 2n/Q,V1- D21024 beschränkt. 
Die Implementierung wird für den Fall einer harmonischen Anregung durch den 
Vergleich der berechneten Lyapunov-Exponenten A; mit den entsprechenden Floquet- 
Multiplikatoren Ar durch den Zusammenhang A; = Re(In(Af)Q,/27) verifiziert. 
Für weitere Informationen zur Theorie und numerischen Berechnung von Floquet- 
Multiplikatoren wird auf [10, 165] verwiesen. 

Um die einzelnen Instabilitätsmechanismen gezielt herauszuarbeiten, werden die 
erhaltenen Ergebnisse zunächst für einen harmonischen und anschließend für einen 
stochastischen Reibwert vorgestellt. Der Fall eines konstanten Reibwerts kann mit einer 
Eigenwertanalyse behandelt werden und führt auf die in Abschnitt 2.3.2 hergeleiteten 
Ergebnisse. Im Falle eines harmonischen Reibwerts us = Acos(Q,T) wird die Stabilitäts- 
karte durch das Auftreten zusätzlicher parametrischer Resonanzen ergänzt. Bei Mehr- 
freiheitsgradsystemen mit Parameteranregung können im Allgemeinen stabilisierende 
und destabilisierende Effekte bei den Fundamental- und Kombinationskreisfrequenzen 
Q./wı = 2/n, O,/w2 = 2/n und Q./(wı + w2) = 1/n fürn = 1,2,... auftreten, die 
von den beiden Eigenkreisfrequenzen om und wz des Systems (5.7) ohne parametrische 
Anregung abhängen. Welche Auswirkungen die Parameteranregung im spezifischen 
Fall hat, hängt dabei stark von der Gleichungsstruktur ab. Aus diesem Grund wird das 
homogene Gleichungssystem (5.7) durch die Modaltransformation (Ax, Az)" = Rq in 
Abhängigkeit der modalen Koordinaten q = (91, q2)" 


2 
"1 


_ Ha -rııfaı Tata, 
te Ach 27 "äech 
c 


RÉI rura 2 
21 
dargestellt. Die Modalmatrix 
Ra i | 
To 72 
enthalt dabei die beiden auf die Massenmatrix M normierten Eigenvektoren des M-K 
Systems aus Gleichung (5.7). Für den Fall gleichmäßiger Dämpfung 6, = 62 = ô 
ergeben sich die Dampfungsbeiwerte 61; = 622 = ô und 6, = 0. Für das vorliegende 


System lässt sich in der Matrix für die Amplituden der Parameteranregung jeweils die 
Phasenverschiebung 7 auf den Haupt- und Nebendiagonalen feststellen. 
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O; A/V 
2201, a 0.25 
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(a) Stabilitatskarte für eine harmonische Anregung (b) Stabilitätskarte fiir eine harmonische Anregung 
mit A = 0.05Y2 und va = Ui < Hi 


Abbildung 5.3: Stabilitätskarten für das homogene Gleichungssystem (5.7) und einen harmonischen Reibwert 
Us = Acos(Q-r) und ô1 = ö2 = 0 - dunkel eingefärbte Flächen kennzeichnen die instabilen 
Bereiche 


Für eine solche Gleichungsstruktur kann aus analytischen Überlegungen gefolgert 
werden [121]: 


e Für Au = ô = 6 existiert ein destabilisierender Effekt bei den Fundamental- und 


Differenzenkombinationskreisfrequenzen.! 


e Fur 6, # 62 lässt sich zusätzlich ein stabilisierender Effekt bei den Summenkombi- 
nationskreisfrequenzen feststellen. Dieser wird auch als Antiresonanz bezeichnet. 


Diese zwei Fälle sollen im Folgenden separat betrachtet werden. Zunächst wird der 
größte Lyapunov-Exponent für den Fall 6; = 62 berechnet, bevor der Fall 61 # 62 
diskutiert wird. 

Für den ungedämpften Fall 5; = 52 = 0 und die Anregungsamplitude A = 0.05¥2 sind 
die berechneten Lyapunov-Exponenten in Abbildung 5.3(a) dargestellt. Dafür wurden 
jeweils die Kreisfrequenz Q, und der deterministische Reibwertanteil ua variiert. 
Dunkel eingefärbte Flächen kennzeichnen einen positiven Lyapunov-Exponenten und 
damit Instabilität. Diese tritt, wie erwartet, bei den Fundamental- und Differenzenkom- 
binationskreisfrequenzen auf. Die bereits aus Abschnitt 2.3.2 bekannte Stabilitätsgrenze 
Ha = ui unterteilt den Parameterraum in die im wesentlichen instabilen und stabilen 
Bereiche ua 2 ui, wobei sich nun aufgrund der Parameteranregung auch kleinere stabile 
Bereiche für ua > u finden lassen. Der Einfluss der Anregungsamplitude A auf die 
instabilen Äste ist in Abbildung 5.3(b) für ua = uı < ui dargestellt. Dafür wurden jeweils 
die Kreisfrequenz Q, und die Standardabweichung R, = A/ V21 + 3 variiert. Es fällt auf, 


1 Fiir diesen Fall lässt sich für das Gleichungssystem (5.7) alternativ auch die Koordinatentransformation 


Ax = qı exp(—6t) und Az = q2 exp(—6T) nutzen. Hierdurch wird die Dämpfungsmatrix eliminiert und 
auf den Hauptdiagonalen der Steifigkeitsmatrix entstehen die zusätzlichen Einträge —67. 
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dass die Instabilitatsbereiche insbesondere bei den doppelten Eigenkreisfrequenzen 
und bei den Differenzenkombinationskreisfrequenzen ausgeprägt sind. Letzteres ist 
eher untypisch fiir mechanische Systeme. 

Um die Auswirkungen des stochastischen Reibwerts für den Fall der gleichmäßigen 
Dämpfung 6, = 62 > 0 zu untersuchen, werden die Lyapunov-Exponenten bei ua = 
UU und 6; = 62 = 0.01 für harmonische und stochastische Parameteranregungen 
mit D = 0,0.01,0.35 berechnet. Die Ergebnisse sind fiir die Standardabweichung 
R, = 0.05 in Abbildung 5.4(a) und für die Standardabweichung R, = 0.1 in Abbil- 
dung 5.4(b) ausgewertet. Instabilität mit Wahrscheinlichkeit 1 tritt für einen positiven 
Lyapunov-Exponenten auf, asymptotische Stabilität mit Wahrscheinlichkeit 1 lässt sich 
für negative Lyapunov-Exponenten festhalten. Die stochastische Parameteranregung 
weicht die Grenzen der ursprünglich klar separierbaren destabilisierten Bereiche 
auf. Dies kann auf das nichtverschwindende Leistungsdichtespektrums (5.6) bei den 
Fundamental- und Differenzenkombinationskreisfrequenzen zurückgeführt werden. 
Während der destabilisierende Effekt bei den doppelten Eigenkreisfrequenzen für 
D = 0.01 noch erhalten bleibt, überdauert für das höhere Dämpfungsmaß D = 0.35 nur 
der destabilisierende Effekt bei den Differenzenkombinationskreisfrequenzen. 

Eine qualitativ gleiche Entwicklung ergibt sich auch für den Fall 6; # 62 in den 
Abbildungen 5.4(c)-(f). Der erwartete stabilisierende Effekt der Parameteranregung 
ist zumindest bei der größten Summenkombinationskreisfrequenz klar erkennbar. 
Die für eine harmonische Anregung noch deutlich abgrenzbaren stabilisierten oder 
destabilisierten Bereiche weichen durch die stochastische Anregung erneut auf, für 
ein hohes Dämpfungsmaß scheint nur noch der destabilisierende Effekt bei den Diffe- 
renzenkombinationskreisfrequenzen zu überdauern. Des Weiteren ist das Ausbleiben 
einer der beiden destabilisierenden Effekte bei den doppelten Eigenkreisfrequenzen 
interessant, das jedoch nur für D = 0.01 stattfindet. 


Die gezeigten Stabilitätskarten offenbaren das Auftreten von Parameterresonanzen auf- 
grund eines stochastischen Reibwerts. Als Ursache für reibungserregte Schwingungen 
scheint dieses Phänomen noch nicht in der Literatur bekannt zu sein. Interessant ist 
insbesondere der für mechanische Systeme eher untypische destabilisierende Effekt bei 
den Differenzenkombinationskreisfrequenzen, der das Auftreten von reibungserregten 
Schwingungen bei niederfrequenten Anteilen des Reibwerts trotz hohen Systemei- 
genkreisfrequenzen erklären könnte. Darüber hinaus sind aber auch die gefundenen 
Antiresonanzen bei Summenkombinationskreisfrequenzen bei ungleichmäßiger Dämp- 
fung eher kontraintuitiv. Die auftretenden Effekte lassen sich auf den harmonischen 
Anregungsfall zurückführen und hängen von der Wahrscheinlichkeit ab, mit der die 
jeweilige Anregungsfrequenz im stochastischen Prozess auftaucht. Inwieweit der hier 
gefundene Instabilitätsmechanismus bei reibungserregten Schwingungen in realen 


120 


5.2 Stochastischer Reibwert 


—D=0 —D =0.01 
—D =0.35 


A1105? 
2 + 


W1-W2W2 


on 2w 


Zon 


—2 
@1-W2 @1+W2 
= — @W1+W2 

Q: 


(a) Größter Lyapunov-Exponent für Ry = 0.05 und 
61 = 62 = 0.01 


A119? 
21 —D=0 —D=0.01 
1 —D = 0.35 


17-0202 


on 2w 


@W1+W2 


Qr 


201 


DCD w1+@2 


(c) Größter Lyapunov-Exponent für R, = 0.05 und 
ö1 = 262 = 0.02 
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(b) Größter Lyapunov-Exponent fiir Ry = 0.1 und 
61 = 62 = 0.01 
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(d) Größter Lyapunov-Exponent für Ry = 0.1 und 
61 = 262 = 0.02 
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(e) Größter Lyapunov-Exponent für Ry = 0.05 und 
261 = 62 = 0.02 


(£) Größter Lyapunov-Exponent für Rq = 0.1 und 
261 = 62 = 0.02 


Abbildung 5.4: Größter Lyapunov-Exponent von dem homogenen Gleichungssystem (5.7) für einen har- 
monischen Reibwert (D = 0) und für einen stochastischen Reibwert mit D = 0.01, 0.35 bei 


Ha = H1 < Mi 
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Anwendungen zu beobachten ist, kann im Rahmen dieser Arbeit jedoch nicht endgiiltig 
beantwortet werden. Da Parameterresonanzen zumindest bei fehlender oder geringer 
äußerer Dämpfung bereits bei sehr kleinen Anregungsamplituden auftreten, sollten 
die Anregungsintensitäten des stochastischen Reibwerts prinzipiell ausreichen. Es 
wird jedoch vermutet, dass andere Instabilitätsquellen in technischen Systemen mit 
Kontakten dominieren. Im Rahmen von Stabilitätsuntersuchungen zu reibungserregten 
Schwingungen aufgrund eines normalverteilten stochastischen Reibwertanteils können 
zudem qualitativ gleiche Phänomene gefunden werden [271]. 
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Im Rahmen dieser Arbeit wird ein thermomechanisches Modell fiir den Kontakt 
zweier gleitender Körper mit rauen Oberflächen entwickelt und für verschiedene Kon- 
taktkonfigurationen ausgewertet. Das Simulationsmodell ermöglicht die Berechnung 
verschiedener kontaktmechanischer Größen wie beispielsweise den Verschiebungen, 
den Spannungen, dem Temperaturfeld und der tatsächlichen Kontaktfläche. Es ist 
außerdem in der Lage, die Einflüsse von temperaturabhängigen Materialparametern 
abzubilden. Die vorliegende Arbeit schließt damit eine in der Literatur noch vorhandene 
Modellierungslücke und kann als Erweiterung der bisher existierenden Kontaktmodelle 
angesehen werden. Weiterführende Untersuchungen hinsichtlich den zugrundeliegen- 
den Instabilitätsmechanismen bei reibungserregten Schwingungen unter Einbezug 
der Erkenntnisse aus den thermomechanischen Simulationen schließen sich in der 
Folge an. Diese zeigen bisher noch nicht in der Literatur bekannte stabilisierende 
und destabilisierende Auswirkungen des stochastischen Reibwerts in Systemen mit 
tribologischen Kontakten, die auf eine stochastische Parameteranregung zurückgeführt 
werden können. 


Für das Modell werden Einflussfunktionen für die voll-gekoppelten und quasi-statischen 
Gleichungen eines thermoelastischen Festkörpers hergeleitet. Dies gelingt durch die 
Anwendung der Neuber-Papkovich-Potentiale und der Verwendung des im Kontext der 
Kontaktmechanik bereits etablierten Halbraummodells. Die auf diese Weise erhaltenen 
analytischen Ausdrücke beschreiben den Zusammenhang zwischen mechanischen 
oder thermischen Lasten innerhalb eines Rechteckelements auf der Oberfläche und 
den Verschiebungen, den Spannungen und der Temperatur an einem beliebigen Punkt 
des Körpers. Mithilfe der erhaltenen Einflussfunktionen können beliebige räumliche 
oder transiente Oberflächenlasten mittels einer geeigneten Superposition simuliert 
werden. Durch die Verwendung der diskreten, zyklischen Faltung im Bildbereich ist 
die Auswertung der räumlichen und zeitlichen Superpositionen speicherplatzsparend 
und rechenzeiteffizient. Eine Verifikation der hergeleiteten Einflussfunktionen und 
der Implementierung erfolgt über den Vergleich mit analytischen Lösungen aus 
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der Literatur. Infolgedessen können Abschätzungen getroffen werden, inwieweit die 
Berücksichtigung des Gough-Joule-Effektes für die Kontaktmechanik von Bedeutung 
ist. Alle hergeleiteten Zusammenhänge behalten ihre Gültigkeit auch bei Anwendungen 
im Rahmen der Poroelastizität, da dort die zugrundeliegenden Gleichungen strukturell 
identisch sind. Eine detaillierte Beschreibung des Gesamtmodells führt schließlich 
zu einem zeitabhängigen, linearen Gleichungssystem mit Ungleichungsnebenbe- 
dingungen, das durch ein CG-Verfahren in jedem inkrementellen Zeitschritt gelöst 
wird. Eine Auswertung des Simulationsmodells erfolgt für isotrope und anisotrope 
Oberflächenkonfigurationen insbesondere in Hinblick auf die Auswirkungen der 
zweiseitigen Rauigkeit. Die verwendeten, künstlich erzeugten Oberflächen besitzen 
dabei sowohl ein vorgegebenes Leistungsdichtespektrum als auch eine vordefinierte 
Verteilungsdichte. Die Auswertung des thermomechanischen Kontaktmodells führt 
letztendlich zu folgenden Erkenntnissen: 


e Die Berechnung des Reibwerts für zwei gleitende Körper mit rauen Oberflächen 
bei unterschiedlichen Relativgeschwindigkeiten und Lasten zeigt, dass der Reib- 
wert von ebenjenen Größen abhängt. Die Geschwindigkeitsabhängigkeit wird 
hauptsächlich durch die Temperaturentwicklung und die Versagensscherspan- 
nung beeinflusst, während sich die Schiefe auf die Lastabhängigkeit auswirkt. 
Für äquivalente Oberflächen, die im Mittel normalverteilte Höhenwerte besitzen 
(und somit E[Ssx] = 0), lässt sich ein lastunabhängiger Reibwert auch bei einer 
Relativbewegung der Kontaktkörper feststellen. Da die Verteilungsdichte der 
Oberflächenhöhen durch den Einlaufprozess verändert wird, hat dieser einen 
signifikanten Einfluss auf die Ergebnisse. Die erhaltenen Erkenntnisse zum sto- 
chastischen Reibwertanteil sprechen zudem gegen eine Modellierung desselben 
als weißes Rauschen. 


e Die konsequente Modellierung der Oberflächenrauigkeiten von beiden Kontakt- 
körpern offenbart ein durchweg transientes Kontaktgeschehen. Dazu zählt in 
besonderem Maße die ebenfalls in Experimenten beobachtbaren Blitztemperatu- 
ren, aber auch der stochastische Anteil des Reibwerts sowie die am Reibprozess 
beteiligten Oberflächenbereiche. Modellierungen, die lediglich eine statische 
äquivalente Oberfläche beinhalten, sind nicht in der Lage, diese systemspezi- 
fischen Eigenschaften in gleichem Maße abzubilden. Des Weiteren finden sich die 
isotropen und anisotropen Eigenschaften der Oberflächenstrukturen im Reibwert 
wieder. 


e Die Berücksichtigung von Reibungswärme und thermoelastischen Verformungen 
deutet einerseits auf einen generell vernachlässigbaren Einfluss des Gough-Joule- 
Effektes und andererseits auf einen kleinen Einfluss der thermoelastischen Verfor- 
mungen bezüglich des Reibwerts hin. Darüber hinaus sollte die Temperaturent- 
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wicklung in der Modellierung bei gleitenden Kontakten berticksichtigt werden, da 
diese unmittelbaren Einfluss auf die Materialparameter hat. In Kombination mit 
der zweiseitigen Rauigkeit konnen Blitztemperaturen erstmals richtig modelliert 
und die auftretenden transienten Verschiebungs- und Spannungsfelder korrekt 
dargestellt werden. Der Einbezug von temperaturabhängigen Materialparametern 
ist möglich. 


Im Rahmen mehrerer Abschlussarbeiten wurde die Modellierung um zwei wesentliche 
Aspekte erweitert. Zum Ersten vermeidet die Verwendung von nicht-äquidistanten 
Gittern die größtenteils unnötig feinen örtlichen oder zeitlichen Diskretisierung [277]. 
Die Anwendung auf partiell gleitende Kontakte knüpft anschließend die Verbindung 
zwischen den Zuständen Haften und Gleiten, sodass eine vollständige Parametrisie- 
rung von dynamischen Reibgesetzen in Kombination mit den thermomechanischen 
Kontaktsimulationen erfolgen kann. Die Ergebnisse der Abschlussarbeit deuten jedoch 
darauf hin, dass die Verwendung von nicht-äquidistanten Gittern zumindest bei 
gleitenden Konfigurationen noch keine zwingenden Vorteile mit sich bringt. Dies liegt 
einerseits an höheren Rechenzeiten und andererseits an der zwangsläufig auftretenden 
Problemstellung der dynamischen Vernetzung. 

Zum Zweiten führt die Verbindung des Multiple-Scales-Verfahrens mit der Reynolds- 
gleichung auf die Modellierung von geschmierten Kontakten mit zweiseitiger Rauig- 
keit [274, 275]. Daraus folgen sogenannte Flusskoeffizienten, die aus repräsentativen 
Ausschnitten der Oberflächen erhalten werden und die eine Abbildung der durch 
die Rauigkeit verursachten Effekte in die makroskopische Skala ermöglichen. Durch 
diese Vorgehensweise lassen sich Stribeck-Kurven für eine Vielzahl unterschiedlicher 
Kontaktkonfigurationen berechnen. Während die zweiseitigen Oberflächenrauigkeiten 
im Bereich der Flüssigkeitsreibung einen geringen Einfluss auf die Stribeck-Kurve 
aufweisen, sind die komplexen Wechselwirkungen im Bereich der Mischreibung 
noch nicht vollständig verstanden. Hier können weitere physikalische Effekte wie 
beispielsweise Kavitation oder Temperaturentwicklung mit einbezogen werden. 


Untersuchungen zu den Auswirkungen der gefundenen Reibwerteigenschaften im Kon- 
text von reibungserregten Schwingungen schließen sich den Auswertungen des thermo- 
mechanischen Simulationsmodells an. Die Untersuchungen konzentrieren sich jeweils 
auf die Stabilität eines stationären Systemzustands und greifen auf ein Minimalmodell 
mit zwei Freiheitsgraden zurück. Die separate Betrachtung von deterministischen und 
stochastischen Anteilen des Reibwerts und die Wahl jeweils möglichst allgemeiner 
Parametrisierungen ermöglichen es, die zugrundeliegenden Instabilitätsmechanismen 
gezielt herauszuarbeiten. Die vorgenommenen Modellierungen und Simulationen 
führen zu folgenden Erkenntnissen: 
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e Die physikalisch sinnvolle Modellierung und Parametrisierung des stochastischen 
Reibwertanteils gelingt durch nichtlineare Filtergleichungen. Der erhaltene sto- 
chastische Prozess genügt der vorgegebenen Autokorrelationsfunktion als auch 
der Verteilungsdichte. Eine Modellierung und Simulation des stochastischen 
Reibwertanteils in beliebigen dynamischen Systemen ist somit durch das Anhän- 
gen der Filtergleichungen problemlos möglich. Es resultieren jedoch stochastische 
Differentialgleichungen, die dienumerische Behandlung erschweren. Methodisch 
ist die Vorgehensweise identisch zur Verwendung der bereits weit verbreiteten 
und gut entwickelten dynamischen, aber deterministischen Reibgesetze. Eine 
Modellierung des stochastischen Reibwertanteils als weißes Rauschen sollte 
vermieden werden. 


e Die Untersuchung der Auswirkungen etwaiger stochastischer Anteile des Reib- 
werts im Kontext von reibungserregten Schwingungsphänomenen offenbart in 
diesem Zusammenhang einen Instabilitätsmechanismus, der so noch nicht in der 
Literatur bekannt ist. Die auftretenden, stabilisierenden und destabilisierenden 
Effekte der stochastischen Parameteranregung lassen sich anhand der Gleichungs- 
struktur und dem Leistungsdichtespektrum erklären. Interessant ist insbesondere 
die gefundene Destabilisierung bei Differenzenkombinationskreisfrequenzen, die 
für mechanische Systeme eher untypisch ist. 


Die Untersuchungen wurden im Rahmen einer weiteren Abschlussarbeit auf schiefe 
Anregungsprozesse erweitert [276]. Der Anregungsprozess wird dabei durch die Metho- 
den aus Unterkapitel 2.1 erzeugt. Die Ergebnisse lassen auf einen zwar vergleichsweise 
kleineren, aber dennoch vorhandenen Einfluss der Verteilungsdichte schließen. 


Es bieten sich viele Ansatzpunkte für weiterführende Forschungsarbeiten an. So 
können die für die Kontaktsimulationen künstlich generierten Oberflächenhöhen zwar 
bereits für beliebige Leistungsdichtespektren und Verteilungsdichten rechenzeiteffizi- 
ent simuliert werden, dies ist jedoch aktuell auf periodische Oberflächen beschränkt. 
Entsprechende Filterfunktionen würden eine Oberflächenerzeugung während der 
eigentlichen Kontaktsimulation ermöglichen und wären damitin der Lage, realistischere 
Kontaktsimulationen durchzuführen. Für die auf diese Weise simulierten technischen 
Oberflächen sind sicherlich auch die (zum aktuellen Zeitpunkt noch nicht verfügbaren) 
Ergebnisse der Surface Topography Challenge interessant [112]. 

Die Erweiterung des vorgestellten thermomechanischen Kontaktmodells um plastische 
Verformungen inklusive Verfestigung kann ebenso erfolgen, wie die Berücksichtigung 
von oberflächennahen Schichten und die Modellierung von chemischen Reaktionen. 
Simulationsmodelle für die einzelnen Bereiche sind in der Literatur bereits größ- 
tenteils bekannt, eine sinnvolle Kombination aller Teilmodelle steht allerdings aus. 
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6 Zusammenfassung und Ausblick 


Darüber hinaus wäre die Einarbeitung von weiterführenden Verschleißmodellen, die 
physikalisch fundierte Modellierungen der unter Umständen gleitrichtungsabhängigen 
Versagensscherspannung sowie in diesem Zusammenhang auch die Erweiterung um 
anisotrope Materialgesetze aufschlussreich. 

Für die aktuell gewählte Randbedingung einer konstanten Einrückung (oder einer kon- 
stanten Normalbelastung) bieten sich außerdem realistischere Alternativen an. Denkbar 
ist beispielsweise ein aus Experimenten bekannter Einrückungsverlauf oder auch eine 
skalenübergreifende Simulation mit gleichzeitiger Berücksichtigung von makrosko- 
pischer Dynamik und mikroskopischem Kontakt. Zumindest Letztere benötigt beim 
aktuellen Stand der Forschung jedoch noch zu lange Rechenzeiten. Auch sollten weitere 
Forschungsbemühungen hinsichtlich der Weiterentwicklung numerischer Methoden 
unternommen werden. 

Bezüglich der Modellierung des stochastischen Reibwerts lässt sich der Zusammenhang 
zwischen der Varianz des Rauschprozesses und der Last sowie die Abhängigkeit 
von der tatsächlichen Relativgeschwindigkeit ergänzen. Weiterführende analytische 
Stabilitätsuntersuchungen können erfolgen. Da sich die vorliegende Arbeit zudem auf 
die Stabilität von Ruhelagen beschränkt, ist eine Erweiterung der Untersuchungen auf 
Grenzzyklen in Betracht zu ziehen. Dazu gehört auch der Übergang der Zustände 
Haften und Gleiten oder die Auswirkungen anisotroper Reibgesetze [172, 267, 268]. 
Bereits deterministische Analysen stoßen hierbei auf unzählige dynamische Schwin- 
gungsphänomene. In diesem Kontext wären auch numerische Integrationsverfahren 
für stochastische Differentialgleichungen mit variabler Schrittweite hilfreich. 

In allen Bereichen steht eine experimentelle Validierung aus. In diesem Zusammenhang 
sind systematische Vermessungen und Auswertungenbei tribologischen Systemen auch 
hinsichtlich stochastischer Eigenschaften hervorzuheben. Ein erster Ansatz wäre hier 
der Vergleich von den Ergebnissen des vorliegenden thermomechanischen Kontaktmo- 
dells mit Experimenten bei simultaner Vermessung der beteiligten Oberflächen vor, 
während und nach den Versuchen. 
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Anhang 


A Poroelastizitat 


Ein mit Fluid durchsetztes, poröses Medium kann durch die beiden voll-gekoppelten 
Gleichungen von Biot [29] 


3A + 2u 
wAr+(A+u)V(V-r)= 3H 


BA +2u 9, -(3 STE 


Vo, (A.1) 


kAo - (A.2) 


3H at "Ip sm? Ia 


beschrieben werden. Sie enthalten neben dem Permeabilitatskoeffizienten k und dem 
Flüssigkeitsdruck o die beiden Koeffizienten 1/H und 1/R, die die Kompressibilität 
des Mediums beziehungsweise die Änderung des Flüssigkeitsgehalts messen. Die 
Gleichungen (A.1) und (A.2) stimmen strukturell mit den homogenen Gleichungen für 
den thermoelastischen Festkörper (2.15) und (2.16) überein. Ein vollständig gesättigtes 
Medium kann durch die Substitution von H = R = (3A +2 u) /3 [29] beschrieben werden, 
wodurch sich das Gleichungssystem zu 


uAr +(A+u)V(V-r)=Vo, 


d 
kAo - —V-r= 
o ER r =0 


ändert. Für einen äquivalenten Effekt bei den thermoelastischen Gleichungen ist für 
die Temperaturleitfähigkeit x — oo notwendig, es folgen direkt € — oo und 


A+2u)K 
lim x. = lim ze ( 2 
K>% vele (31+2u)a2Tı 


sowie die Péclet-Zahl 


_ 4a2— 4a? (BA +2u) a?Ty 
lim Pe = lim = 
Ko K>% K-At (A + 2u) KAt 
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B 


Einflussfunktionen 


Für mechanische Oberflächenbelastungen ergeben sich die dimensionslosen Einfluss- 
funktionen aus Abschnitt 3.1.1 zu 


mit 


Cry = Dzy|x=2e41 + Prx|x=22-1 - ®zx|x=2z-1 Olai, 
Y=2n+1 Y=2n-1 Y=2n+1 Y=2n-1 


we 


zy = Dzy|x=241 + Tele: =26-1 — zy|x=2z-1 - Dzy|x=2e41, 
Y=2n+1 =2n-1 Y=2n+1 Y=2n-1 


Ca = Ozz|xn2E41 + Dzz|x=08-1 — Pez |xare-1 — Ozz|x=2e41, 
Y=2n+1 Y=2n-1 Y=2n+1 Y=2n-1 


Tzx = Ozx|x=2¢41 + Or |x=22-1 - Ozx|x=2¢-1 - Oxx|x=2841, 
Y=2n+1 Y=2n-1 Y=2n+1 Y=2n-1 


Tay = Oxy |¥- SE + Ozy|x= = ı- Ozy|x=2¢-1 - ©zy|x= SZ 
Y=2n+1 -1 Y=2n+1 Y=2n-1 


Tzz = ©zz|x=2z+1 + ney - O2 |x-2:-1 - Ozz|x=2641, 
Y=2n+1 Y=2n-1 Y=2n+1 Y=2n-1 


Q, = D,|x-2¢41 + By|x=22-1 - als - Og|x=2e41 
Y=2n+1 Y=2n-1 Y=2n+1 Y=2n-1 


Y 
d =ßY In (x? + py’) + 2X arctan (=| 


d =XIn (x? n PY) + 2BY arctan Di , 


d == In D + 4/X? nl + CG In (x +4/X2 + pae) : 


Q., =In Le +4/X2 + er , 
Oz, =In [-x +4/X2 + er , 
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B Einflussfunktionen 


(S =] sign (X) sign (BY), 


arog 


P, =- gY 


Für thermische Oberflächenbelastungen und vorgegebenem Wärmestrom ergeben 
sich die dimensionslosen Einflussfunktionen aus Abschnitt 3.1.2 


T, = =D, 3 Gel Säi 7 Oo X=2é-1 — Oe X=2é-1 — Oo vel dt 
Y=2n+1 Y=2n-1 Y=2n+1 Y=2n-1 
Ne 
mit 
= X BY 
Os = erf Le erf Li A 
Die Substitution t = 4/5 I wird genutzt um die Singularität zu umgehen. Für die 


diskreten Zeitpunkt t, = hAt und t = EA ändern sich die Integralgrenzen zu Vk - h 
und Vk. 
Außerdem ergibt sich 


ae 
AE 


Sa 4 
zz = P. 08zz|x=2£+1 + 08zz|x=2£-1 = 08zz|x=2£-1 — O08zz|x=2&+1 dt 
S E Y=2n-1 Y=2n+1 Y=2n-1 
h 
At 


mit 


Ozz = — XBY 


d baie a, L erf al ad Pe. 


Die Substitution 7 = = wird genutzt. Für die diskreten Zeitpunkte tp = hAt und 
t = kAt werden die Integrationsgrenzen zu k — h und k. 
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B Einflussfunktionen 


Für thermische Oberflächenbelastungen und vorgegebenem Temperaturfeld ergeben 
sich die dimensionslosen Einflussfunktionen aus Abschnitt 3.1.2 


Ts = Og|x=2e+1 + Og|x=2e-1 - Og|x=2¢-1 - Oslx=2cH1 
Y=2n+1 Y=2n-1 Y=2n+1 Y=2n-1 


mit 
T ; 
O5 = z Sign (X)sign (BY). 


Außerdem ergeben sich 


t 
8 
=— O@zx|x=2841 + Oozx|x=2¢-1 — C@zx|x=2-1 — Oozx|x=2¢+1] AT, 
Y=2n+1 Y=2n-1 Y=2n+1 Y=2n-1 


AF 
cÉ 
Ar 

f Lech =2641 + Ogzy|X=2E-1 — O@zy|X= oe — Oezy [x= scha 


=2n+1 Y=2n-1 Y= Y=2n-1 


mit 


ei 
sl 
H 
ma 


Bet _ X? + pY? 
wen (VP zl KEENT dee Dr | 
e x = Pe) 


“Se Ter — 


1 X 
ozy = erf | VPe—Je P, m 
zy pY? | Er) e X2 +4 BY? 


Die Substitution t = wird genutzt. Für die diskreten Zeitpunkte tp = hAt und 


t = kAt werden die Integrationsgrenzen zu k — h und k. 
Des Weiteren ergibt sich 


a 


— 2 
WS welx=25+1 + welx=25-1 - welx=28-1 - We|xX=2E41] dT 
e Y=2n+1 Y=2n-1 Y=2n+1 Y=2n-1 
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B Einflussfunktionen 


mit 
G [X2 + BY? e je: [X2 + BY? i 
We Sg er en: T 


P Das BY\ P x 
erf BI e “16? ert (vP) e "ie 


Die Substitution t = H wird genutzt. Für die diskreten Zeitpunkte t, = hAt und 
t = kAt werden die Integrationsgrenzen zu Vk — h und Vk. 
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C  Fundamentallésung der 
Wärmeleitungsgleichung 


Die Wärmeleitungsgleichung 


00 
AO Te yt =O )= 
K E ost UI 0 


wird für den instantanen Wärmeeintrag 4 zum Zeitpunkt t = 0 im Ursprung des 
Vollraums gelöst. Aufgrund der Symmetrie bietet sich zunächst die Auswertung der 
Wärmeleitungsgleichung in Kugelkoordinaten 


Ke 0 [,,00\ 90 =, Pees NER 
rl E) Ve Fa 


an, bevor die Laplacetransformation 


+00 


L{O(R, t)}(R,s) = d O(R, edt = 6 


0 


unter Berticksichtigung der Anfangsbedingung die Differentialgleichung 


Ke d'Long ` a 
rl =) d 


liefert. Deren allgemeine Lösung 


g = Ne EVE LZAVE, 
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C Fundamentallösung der Wärmeleitungsgleichung 


kann unter Berücksichtigung der Nebenbedingung Jim 6* = 0 und erfolgter Rücktrans- 


formation [193] als 


geschrieben werden. Die noch verbliebene Konstante cı lässt sich aus der Bilanz der zu- 
oder abgeführten Wärmemenge durch das Volumenintegral 


K 
— | 6dv=4q 
Ke RER 


zu cı = 4/4rıK bestimmen, sodass die Fundamentallösung 


= x2+y2+z2 


Kee Jet 


8K (1K et)?” 


Č = 


erhalten wird. Sie besitzt für den Zeitpunkt des Wärmeeintrags eine Singularität im 
Ursprung. 
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D Stochastische Stabilitätstheorie 


Ausgangspunkt sei eine beliebige stochastische Differentialgleichung im Sinne von Itö 
dX,=a(X,,s)ds+b(X,,s)dW,, (D.1) 


wobei sich die nachfolgenden Ausführungen an Khasminskii [123], Arnold [11] sowie 
Lin und Cai [139] orientieren. Die Stabilitätsaussagen beziehen sich auf die triviale 
Lösung X, = 0, dementsprechend muss für den Driftterm a (0,s) = 0 und für den 
Diffusionsterm b (0, s) = 0 gelten. 

Die Lösung X; = 0 heißt stabil mit Wahrscheinlichkeit 1 für s > so, wenn für beliebige 
so >Ounde >0 


li X = 
Jim, P {sup | sll > d 0 


EZE 


gilt. Des Weiteren heift die Lösung X, = 0 asymptotisch stabil mit Wahrscheinlichkeit 
1, wenn sie stabil mit Wahrscheinlichkeit 1 ist und zusätzlich 


lim P {lim IX: = o} =1 
X so —0 so 
gilt. Ist diese Eigenschaft unabhängig von den Anfangswerten 
P | lim Xl] = ol =1, 
Ek deal 


wird sie als asymptotisch stabil mit Wahrscheinlichkeit 1 im Ganzen bezeichnet. 
Die Definitionen gehen für einen verschwindenden Diffusionsterm b (X;,s) = 0 in 
diejenigen des deterministischen Falls über, sodass die stochastische Stabilitätsdefinition 
als Verallgemeinerung angesehen werden kann. 

Die Stabilitätsaussagen können für die triviale Lösung X; = 0 des Systems (D.1) 
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D Stochastische Stabilitätstheorie 


aufgrund des Hartmann-Grobmann Theorems für hyperbolische Fixpunkte mithilfe 
des größten Lyapunov-Exponenten 


1 AX 
Aı = lim — In WAX 
ees |[AXs9| 


und den linearisierten Systemgleichungen 


_ da(x,s) AX.ds + ob (x, 8) 
dr dr 


x=0 


AX,dW; 
x=0 
getroffen werden (Theorem 7.1 und Theorem 7.2 von Khasminskii [123]). Die Existenz 
des Grenzwertes ist durch das Multiplikative Ergodentheorem von Oseledec gesichert. 
Die Lösung X, = 0 kann fiir Au < 0 als asymptotisch stabil mit Wahrscheinlichkeit 
1 und ftir A; > 0 als instabil mit Wahrscheinlichkeit 1 klassifiziert werden. Fur den 
Sonderfall A; = 0 lässt sich auch für lineare Systeme zumeist nur eine asymptotische 
Stabilität ausschließen. 
Neben der hier vorgestellten Variante der stochastischen Stabilitätstheorie existiert 
ebenfalls der historisch ältere Ansatz, die Stabilitätsaussagen an die Erwartungswerte zu 
koppeln. Diese Vorgehensweise lohnt sich in der Regel jedoch nur bei normalverteilten 
Prozessen, da diese komplett über die ersten beiden Momente beschrieben sind. 
Weiterführende Informationen hierzu lassen sich in [11, 123] finden. 


Die Ausdrücke für die Driftterme in Gleichung (5.9) lauten 


f= k -u sin(2p) — 2ô1 oner! cos(9)? - (2 sin(2w) + 2ô2 siny?) sin(9)? 


- (2 Costa) sin(y) + ; sin(p) cos() (C12 — ua — Asin(@) SEIN sin(29), 
fo =(1 — C11) cos(p)? - 1 - 6; sin(2¢) 

— cos(p) cos() tan(9) (C12 - ua — Asin(b) cos(0)) , 
fa =- Cz cos(@) cos(y) cot(S) — C22 cos(W)* — 1 — 62 sin(2), 
faz ( ou = : sin(2p) + 61 sin(p)? - = sin(2W) — 59 sin?) sin(29) 

+ sin(p) cos(w) sin(9)? (C12 — ua — Asin(#) cos(0)) 


-Čz cos(P) sin(y) cos(9)?) sign (sin(2$)), 


f=- SE (2 tan(@) sin(0)? — cot(#) c0s(0)?) sign (sin(2¢)) , 
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D Stochastische Stabilitätstheorie 


2 
14+21 


fs =Q - DQ, sin(20) I + cool. 


Die Ausdriicke fiir die Diffusionsterme in Gleichung (5.9) lauten 


h=- SPor sin(0) sign (sin(ġ)), 


f=- oe cot(&) cos(@) sign (cos(¢)) . 
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Abhängigkeit verschiedener physikalischer Pa- 
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